Smoking and COVID Living Review

0.1 Data preparation

0.1.1 Load packages

• Load packages: load packages and useful scripts for the analysis

library('tidyverse', quietly = T)
library('ggplot2', quietly = T)
library('dplyr', quietly = T)
library('readxl', quietly = T)
library('linelist', quietly = T)
library('googlesheets4', quietly = T)
library('naniar', quietly = T)
library('gridExtra', quietly = T)
library('snakecase', quietly = T)
library('epitools', quietly = T)
library('meta', quietly = T)
library('knitr', quietly = T)
library('kableExtra', quietly = T)
library('snakecase', quietly = T)
library('here', quietly = T)
library('flextable', quietly = T)
library('DT', quietly = T)

0.1.2 Download from Google sheets

• Load data: import raw dataset directly from Google sheets We import the data from the dataset that is available in the “data for analysis” folder on the Google drive. This is available to be viewed with the link and is only available for viewing not editing.

In future versions we may open a tab in the Google sheet to allow people to submit relevant papers for us to incorporate in our review.

sheets_id <- as_sheets_id("https://docs.google.com/spreadsheets/d/1fuOcaN7K_f3msmt0aQAWwOWjgD4SsEt7iJ_0to1wotA/edit?usp=sharing")
search_details <- read_sheet(sheets_id, range = 'article_screening')
data_study_general <- read_sheet(sheets_id, range =  'general_details')
table_1 <- read_sheet(sheets_id, sheet = 'pop_descriptives')
table_2 <-  read_sheet(sheets_id, sheet = 'testing')
table_3 <-  read_sheet(sheets_id, sheet = 'hospitalisation')
table_4 <-  read_sheet(sheets_id, sheet = 'severity')
table_5 <-  read_sheet(sheets_id, sheet = 'mortality')
table_6 <-  read_sheet(sheets_id, sheet = 'updated_quality_appraisal')

The data is downloaded from the sheet and split into dataframes to represent the current sheets and analysis.

0.1.3 Raw Data

The first sheet is Search details we use this sheet to keep track of the number of studies found in our search and keep track of the number of titles and abstracts screened and number of papers included in each review version. It presents the data in groups of search periods.

flextable(head(search_details[,1:11]))

date_screening in format (2020-05-22)	ovid_number_results	ovid_screened	ovid_excluded	ovid_included	medrxiv_number_results	medrxiv_screened	medrxiv_included	other_source_included	total_included	review_version
2020-04-21	15	NA	NA	4	NA	NA	NA	24	28	v1
2020-05-05	57	7	4	3	37	10	0	0	3	v2
2020-05-16	67	12	12	0	70	64	22	6	28	v3
2020-06-02	297	48	44	4	158	158	32	2	38	v4
2020-06-23	477	72	63	6	157	157	41	2	49	v5
2020-07-14	656	30	26	4	105	105	23	1	28	v6

The second sheet is General study details this includes the lead author, where the study is published, the date it was published, the country in which the data originated and the DoI or link to the publication for our records.

This is a searchable list of all the studies that have met out inclusion and exclusion criteria up until -Inf

The third sheet produces Table 1 this contains the details of the population of the studies in our review. Primarily includes smoking status but also the proportion of females in the study and the age of participants.

The fourth sheet produces Table 2 this contains the data extracted from relevant studies on testing data and smoking status. Most studies will not have data on this outcome however we filter these subsequently. These are sparse datasets as many studies do not report on all levels of smoking exposure.

The fourth sheet produces Table 3 this contains the data extracted from relevant studies on hospitalisation and smoking status.

The fourth sheet produces Table 4 this contains the data extracted from relevant studies disease severity data and smoking status. Disease severity is differentially reported across the studies in some there is reference to WHO severity score in others we have based the distinction on clinical characteristics or whether patients were admitted to Intensive or Critical care areas.

The fifth sheet produces Table 5 this contains the data extracted from relevant studies on COVID-19 related mortality and smoking status. A caution needed when interpreting this data is that most studies have not followed patients up for sufficient amounts of time to be able to accurately report on disease associated mortality.

The final sheet was used to collate data on Study quality however, we have discontinued the use of this tool from version 3 due to the large number of components that were unable to be assessed in the rapidly expanding area of COVID-19 related research. We have moved onto a different method of assessing study quality which is explained in greater detail in the manuscript.

0.1.4 Clean data and save updated versions as RDS files

The data from the google sheets document is then cleaned using a script from the linelist package which has been developed by T. Jombat and Z. Kamvar within the R Epidemics Consortium

The script clean_data() ensures consistency among column names, country names, dates and missing values in the included sheets. Data is also saved as .rds files in the /data_clean folder of the project repository.

search_details <- search_details %>%
  clean_data() %>%
  rename(date_screened = date_screening_in_format_2020_05_22) %>%
  write_rds(here::here('data_clean', 'search_details.rds')) ## The title of the date screened column is to communicate among the review team how we need date formatted

protect_columns_1 <- names(data_study_general) %in% 'doi' ## As the DoI's are links we want to preserve their formatting when running the cleaning script

data_study_general <- data_study_general %>%
  rename('notes' = 11, date_published = `date_published_format(2020-05-12)`) %>%
  clean_data(protect = protect_columns_1) %>%
  mutate(study_id = 1:length(lead_author)) %>% ## We create a unique study ID to aid checking of the outputs later
  write_rds(here::here('data_clean', 'data_study_general.rds'))

review_details <- data_study_general %>%
  select(lead_author, date_published, country, review_version, study_id) %>% ## The version of the review in which included studies were incorporated is noted
  write_rds(here::here('data_clean', 'review_details.rds'))

table_1 <- table_1 %>%
  clean_data() %>%
  mutate(., lower_range = sub('\\_.*', '', .$range)) %>%
  mutate(., upper_range = sub('.*_', '',.$range,)) %>% ## We break apart the age range into an upper and lower level
  select(lead_author, sample_size, female_sex_percent, data_source, median_age, mean_age, iqr_lower, iqr_upper, standard_deviation, lower_range,
         upper_range, current_smoker, former_smoker, current_former_smoker, never_smoker, never_smoker_unknown, current_vaper, former_vaper,
         current_former_vaper, not_stated, missing, total) %>%
  filter(!is.na(lead_author)) %>%
  left_join(., review_details, by = 'lead_author') %>%
  write_rds(here::here('data_clean', 'table_1.rds'))


##For the following tables we're only interested in obtaining data for trials which report on the specific outcome of interest

table_2 <-  table_2 %>%
  clean_data() %>%
  filter(data_on_testing == TRUE) %>%
  left_join(., review_details, by = 'lead_author') %>%
  write_rds(here::here('data_clean', 'table_2.rds'))

table_3 <-  table_3 %>%
  clean_data() %>%
  filter(data_on_hospitalisation == TRUE) %>%
  left_join(., review_details, by = 'lead_author') %>%
  write_rds(here::here('data_clean', 'table_3.rds'))

table_4 <-  table_4 %>%
  clean_data() %>%
  filter(data_disease_severity == TRUE) %>%
  left_join(., review_details, by = 'lead_author') %>%
  write_rds(here::here('data_clean', 'table_4.rds'))

table_5 <- table_5 %>%
  clean_data() %>%
  filter(data_on_deaths == TRUE) %>%
  left_join(., review_details, by = 'lead_author') %>%
  write_rds(here::here('data_clean', 'table_5.rds'))


protect_columns_2 <- !names(table_6) %in% 'lead_author'

table_6 <-  table_6 %>%
  clean_data(protect = protect_columns_2) %>%
  left_join(., table_1 %>%
              select(lead_author,
                     not_stated,
                     missing,
                     total),
                          by = 'lead_author') %>%
  mutate(missingness = rowSums(.[8:9], na.rm = T),
         missingness_percentage = (missingness/total)*100) %>% ## Later we want to assess the amount of missingness for each study as it would impact the inclusion of a study into the relevant meta-analyses
  select(-c(not_stated, missing, total, missingness)) %>%
  write_rds(here::here('data_clean', 'table_6.rds'))

a <- data_study_general %>%
  select(lead_author, date_published, source, study_id) %>%
  rename('Lead Author' = lead_author,
         'Date Published' = date_published,
         'Publication Source' = source,
         'Study ID' = study_id)

a$`Lead Author` <- to_upper_camel_case(a$`Lead Author`, sep_out = ", ")
a$`Publication Source` <- to_title_case(a$`Publication Source`)
a$`Publication Source` <- if_else(str_length(a$`Publication Source`) < 5,
                                  toupper(a$`Publication Source`),
                                  to_title_case(a$`Publication Source`))

0.2 Preparing the descriptive analysis

0.2.1 Choosing which data we want to include in the current review?

date_of_update <- Sys.Date()

prev_versions <- c('v1', 'v2', 'v3', 'v4', 'v5')
#This can categorise which studies we want to look at
current_version <- c('v6')
#This will categorise which studies we are including in the current report

analysed_versions <- c('v1', 'v2', 'v3', 'v4', 'v5', 'v6')
#This will incorporate all studies into the current version of the report

We define the current version of the report here. At the moment this needs to be changed directly when compiling the code. This report has been generated for version v6

exclude_from_analysis <- c('isaric_1', 'isaric_2', 'isaric_3', 'miyara_old', 'isaric_4', 'mehra', 'hopkinson', 'gaibazzi', 'miyara_updated', 'russell', "petrilli_old", "magagnoli_old", "niedzwiedz_old", "bello_chavolla_old", "kimmig_old", "russell_old", "zuo_zuo_old", "patel_old")

We don’t want to remove studies from the dataset despite them no longer being included in the analysis. The excluded studies are defined in this chunk of code. We have excluded this version the following studies: isaric_1, isaric_2, isaric_3, miyara_old, isaric_4, mehra, hopkinson, gaibazzi, miyara_updated, russell, petrilli_old, magagnoli_old, niedzwiedz_old, bello_chavolla_old, kimmig_old, russell_old, zuo_zuo_old, and patel_old.

0.2.2 Filter the included studies

Based on the version and excluded studies we produce the tables for analysis. This is a searchable table.

data_study_general <- data_study_general %>%
  filter(review_version %in% analysed_versions) %>%
  filter(!(lead_author %in% exclude_from_analysis)) %>%
  write_csv(., here::here('data_clean', 'data_study_general.csv'))

table_1 <- table_1 %>%
  filter(review_version %in% analysed_versions) %>%
  select(-review_version) %>%
  filter(!(lead_author %in% exclude_from_analysis))

table_2 <- table_2 %>%
  filter(review_version %in% analysed_versions) %>%
  select(-review_version) %>%
  filter(!(lead_author %in% exclude_from_analysis))

table_3 <- table_3 %>%
  filter(review_version %in% analysed_versions) %>%
  select(-review_version) %>%
  filter(!(lead_author %in% exclude_from_analysis))

table_4 <- table_4 %>%
  filter(review_version %in% analysed_versions) %>%
  select(-review_version) %>%
  filter(!(lead_author %in% exclude_from_analysis))

table_5 <- table_5 %>%
  filter(review_version %in% analysed_versions) %>%
  select(-review_version) %>%
  filter(!(lead_author %in% exclude_from_analysis))

a <- data_study_general %>%
  select(lead_author, date_published, source, study_id) %>%
  rename('Lead Author' = lead_author,
         'Date Published' = date_published,
         'Publication Source' = source,
         'Study ID' = study_id)

a$`Lead Author` <- to_upper_camel_case(a$`Lead Author`, sep_out = ", ")
a$`Publication Source` <- to_title_case(a$`Publication Source`)
a$`Publication Source` <- if_else(str_length(a$`Publication Source`) < 5,
                                  toupper(a$`Publication Source`),
                                  to_title_case(a$`Publication Source`))

datatable(a)

a <- flextable(a)%>%
  set_table_properties(width = 1, layout = 'autofit') %>%
  set_caption(caption = 'Studies included in the current analysis')

save_as_docx(a, path = here('data_clean', 'included_studies.docx'))

0.3 Descriptive analysis

0.3.1 The PRISMA diagram will be updated with these numbers

The PRISMA reporting guidelines for systematic reviews the guidance is available here. We use these values to construct the flow diagram of the results from our search strategy and versioning. The numbers included in the quantitative analysis are those meeting quality eligibility criteria and are calculated later in the report.

#Number screened to add to PRISMA
from_prev_version <- data_study_general %>%
  filter(review_version %in% prev_versions) %>%
  filter(!(lead_author %in% exclude_from_analysis)) %>%
  tally()

source(here::here('scripts', 'prisma_function.R'))

PRISMA_v5 <- PRISMA(search_details)
PRISMA_v5 <- as.data.frame(PRISMA_v5)
a <- PRISMA_v5 %>%
  t() %>%
  as.data.frame() %>%
  add_rownames() %>%
  flextable() %>%
  set_header_labels(rowname = 'PRISMA category',
                    V1 = 'Number of studies') %>%
  set_table_properties(width = 1, layout = 'autofit') %>%
  set_caption(caption = 'Data for PRISMA')

save_as_docx(a, path = here('data_clean', 'PRISMA_details.docx'))
a

Data for PRISMA

PRISMA category	Number of studies
ovid_number_results	656
medrxiv_number_results	105
other_source_included	1
titles_screened	762
excluded_abstract	626
full_texts_reviewed	136
full_texts_excluded	108
full_texts_included_narrative	28
total_all_versions	168

0.3.2 We summarise where the studies were conducted

As the pandemic continues studies are being reported from an increasing numbers of countries.

library('ggmap')


#Countries
country <- table_1$country %>%
  to_upper_camel_case() %>%
  recode(., 'Usa' = 'USA', 'Uk' = 'UK', 'SaudiArabia' = 'Saudi Arabia', "SouthKorea" = "South Korea") %>%
  table() %>%
  as.data.frame() %>%
  arrange(desc(Freq)) %>%
  rename('Country' = 1,
         'Number' = 2) %>%
  mutate(Country = as.character(Country))

ggplot(country, aes(x = reorder(Country, desc(-Number)), y = Number))+
  geom_col()+
  coord_flip()+
  theme_bw()+
  labs(title = 'Countries where studies were performed',
       y = 'Number of studies',
       x = 'Country')

a <- flextable(country) %>%
  set_table_properties(width = 0.5, layout = 'autofit')
a

Country	Number
China	45
USA	44
UK	18
Mexico	11
Spain	9
France	8
Italy	6
Multiple	6
Brazil	3
Israel	3
Finland	2
Iran	2
Turkey	2
Australia	1
Canada	1
Chile	1
Denmark	1
Germany	1
Hungary	1
India	1
Kuwait	1
Poland	1
Portugal	1
Saudi Arabia	1
South Korea	1
Sweden	1
Switzerland	1
Thailand	1

country <- subset(country, Country != "Multiple")

country <- mutate_geocode(country, Country)

country_bounds <- country %>%
  select(lon, lat) %>%
  drop_na %>%
  sp::SpatialPointsDataFrame(coords = .[,1:2], data = .) %>%
  sp::bbox() * 1.1

map_world <- map_data(map = "world") %>%
  filter(region != "Antarctica")

save_as_docx(a, path = here('data_clean', 'Countries_performing_studies.docx'))  

0.3.3 A map is produced for where the studies are performed

library(rgdal)
library(leaflet)
library(coronavirus)
update_dataset()

## Updates are available on the coronavirus Dev version, do you want to update? n/Y

my_spdf <- readOGR( 
  dsn= paste0(here::here("data_clean", "world_shape_file")), 
  layer="TM_WORLD_BORDERS_SIMPL-0.3",
  verbose=FALSE
)

from_countries <- c("USA", "UK", "Iran", "South Korea")
to_countries <- c("United States", "United Kingdom", "Iran (Islamic Republic of)", "Korea, Republic of")
country_map <- country %>%
  mutate(Country = plyr::mapvalues(Country, from = from_countries,
                to = to_countries))

from_countries <- c("US", "Iran", "Korea, South", "Czech Republic")
to_countries <- c("United States", "Iran (Islamic Republic of)", "Korea, Republic of", "Czechia")

data("coronavirus")
summary_df <- coronavirus %>% 
  filter(type == "confirmed") %>%
  group_by(country) %>%
  summarise(total_cases = sum(cases)) %>%
  arrange(-total_cases) %>%
  rename("NAME" = country) %>%
  mutate(NAME = plyr::mapvalues(NAME, from = from_countries,
                to = to_countries))

my_spdf@data <- left_join(my_spdf@data, country_map %>%
            rename(NAME = Country),
          by = "NAME") %>%
  left_join(., summary_df %>%
              mutate(total_cases = total_cases/1000),
            by = "NAME") %>%
  mutate(NAME = as.factor(NAME))

my_spdf@data$Number[ which(my_spdf@data$Number == 0)] = NA
my_spdf@data$total_cases[ which(is.na(my_spdf@data$total_cases))] = 0


my_bins <- c(0, 1, 5, 10, 20, 30, 40, 50)
mypalette <- colorBin(palette="YlOrBr", domain=my_spdf@data$Number, na.color="transparent", bins=my_bins)

mytext <- paste(
  "Country: ", my_spdf@data$NAME,"<br/>", 
  "Number of studies included: ", my_spdf@data$Number, "<br/>", 
  "COVID-19 reported cases (thousands): ", round(my_spdf@data$total_cases, 2), 
  sep="") %>%
  lapply(htmltools::HTML)

html_map <- leaflet(my_spdf) %>% 
  addTiles()  %>% 
  setView(lat=10, lng=0 , zoom=2) %>%
  addPolygons( 
    fillColor = ~mypalette(Number), 
    stroke=TRUE, 
    fillOpacity = 0.9, 
    color="white", 
    weight=0.3,
    label = mytext,
    labelOptions = labelOptions( 
      style = list("font-weight" = "normal", padding = "3px 8px"), 
      textsize = "13px", 
      direction = "auto"
    )
  ) %>%
  addLegend( pal=mypalette, values=~Number, opacity=0.9, title = "Number of studies included", position = "bottomright" )

html_map

This map displays the number of studies produced in each country. The number of confirmed COVID-19 cases are shown when resting over each country for reference, data on confirmed cases have been obtained from Rami Krispin’s coronavirus package. The number of cases is updated when the report is generated so will not remain in date, it is current as of 2020-07-20 11:54:57

0.3.4 We sumarise the settings in which they were conducted

Most studies are performed in hospital settings. Some studies are conducted in both such as those in emergency departments where patients are discharged home and followed up. Some registry studies utilising primary care networks are solely conducted in the community.

#Setting
setting <- to_upper_camel_case(data_study_general$study_setting, sep_out = ' ') %>%
  table(.) %>%
  sort('Number of studies', decreasing = T) %>%
  as.data.frame() %>%
  arrange(desc(Freq)) %>%
  rename('Setting' = 1,
         'Number' = 2)

a <- flextable(setting) %>%
  set_table_properties(width = 0.5, layout = 'autofit')
a

Setting	Number
Hospital	119
Community And Hospital	46
Community	7
Not Stated	1
Quarantine Centre	1

save_as_docx(a, path = here('data_clean', 'Countries_performing_studies.docx'))  

0.3.5 We summarise the design of the included studies

The majority were retrospective cohorts. There were fewer case-control and cross-sectional studies. There were some prospective cohorts established and reported.

design <- to_upper_camel_case(data_study_general$study_design, sep_out = ' ') %>%
  recode("Rct" = "RCT") %>%
  table(.) %>%
  sort('Number of studies', decreasing = T) %>%
  as.data.frame() %>%
  arrange(desc(Freq)) %>%
  rename('Design' = 1,
         'Number' = 2)

a <- flextable(design) %>%
  set_table_properties(width = 0.5, layout = 'autofit')
a

Design	Number
Retrospective Cohort	118
Prospective Cohort	34
Cross Sectional	14
Case Control	3
Matched Case Control	3
RCT	2

save_as_docx(a, path = here('data_clean', 'Study_designs.docx'))  

##For supplementary tables
study_design <- data_study_general %>%
  select(lead_author, study_design) %>%
  mutate(study_design = to_upper_camel_case(study_design, sep_out = " "),
         study_design = ifelse(study_design == "Rct", "Randomised Controlled Trial", study_design))

0.3.6 We update the number of participants we have included

• Summary statistics: this table will give us the mean, median and IQR of the samples in the review

#Numbers
a <- summary(table_1$sample_size)
a

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##        6      124      396   116792     1573 17278392

• Total number of participants: 20321806

0.3.7 We want to know if they use a CRF or electronic health record

The majority of data on smoking outcomes is extracted from electronic health records that are routinely collected in the provision of healthcare. Some of the studies utilise specifically designed Case Report Forms (CRFs). The source of information on smoking status is not always reported.

#Source of data
table_1$data_source %>%
  recode('electronic_health_records' = 'Electronic health record',
         'case_report_form' = 'Case report form',
         'not_stated' = 'Not stated') %>%
  table() %>%
  as.data.frame() %>%
  arrange(desc(Freq)) %>%
  rename('Data source' = 1,
         'Number' = 2) %>%
  flextable() %>%
  set_table_properties(width = 0.5, layout = 'autofit')

Data source	Number
Electronic health record	104
Case report form	45
Not stated	25

0.3.8 Now about the smoking

#Studies reporting smokers
current_smok <- table_1 %>%
  filter(., current_smoker != 'NA')

• Current smokers: 98 studies report on the numbers of current smokers

#Studies reporting former smokers
former_smok <- table_1 %>%
  filter(., former_smoker != 'NA')

• Former smokers: 48 studies report on the numbers of former smokers

#Studies reporting never smokers
never_smok <- table_1 %>%
  filter(., never_smoker != 'NA')

• Never smokers: 53 studies report on the numbers of never smokers

#Studies reporting current/former smokers
current_former_smok <- table_1 %>%
  filter(., current_former_smoker != 'NA')

• Current/former smokers: 75 studies report on the numbers of current/former smokers without distinguishing between them, often termed as yes to a smoking history question.

#Studies reporting missing data
missing_smok <- table_1 %>%
  filter(., missing != 'NA') %>%
  mutate(percent = missing/total*100)
minimum <- min(missing_smok$percent)
maximum <- max(missing_smok$percent)
#Studies reporting never/unknown
never_smok_unknown <- table_1 %>%
  filter(., never_smoker_unknown != 'NA')

#Studies with not stated
smok_not_stated <- table_1 %>%
  filter(., not_stated != 'NA')

• Missing, never/unknown and not stated: studies report on the numbers of missing (55), never/unknown (13) or not stated (108) smoking history. This is often framed as absence of a smoking history.
  • Minimum missing: 0.08%
    
  • Maximum missing: 96.36%

#Studies reporting current, former and never smoking status
full_smoking_status <- table_1 %>%
  filter(lead_author %in% current_smok$lead_author) %>%
  filter(lead_author %in% former_smok$lead_author) %>%
  filter(lead_author %in% never_smok$lead_author)

b <- full_smoking_status %>%
  select(lead_author, sample_size, current_smoker, former_smoker, never_smoker, study_id) %>%
  rename('Lead author' = lead_author,
         'Sample size' = sample_size,
         'Current smokers' = current_smoker,
         'Former smokers' = former_smoker,
         'Never smokers' = never_smoker,
         'Study ID' = study_id)

b$`Lead author` <- to_upper_camel_case(b$`Lead author`, sep_out = ", ")

datatable(b,
          caption = "Studies providing data on current, former and never smokers") %>%
  formatRound(columns = c(2:5), digits = 0, interval = 3, mark = ",")

#These are the studies with data on current, former and never smokers

• Studies with complete smoking status: 39

#Studies reporting current or current/former and never smoking
semi_full_smoking_status <- table_1 %>%
  filter(lead_author %in% current_former_smok$lead_author) %>%
  filter(lead_author %in% never_smok$lead_author) %>%
  filter(!lead_author %in% full_smoking_status$lead_author)

b <- semi_full_smoking_status %>%
  select(lead_author, sample_size, current_former_smoker, never_smoker, study_id) %>%
  rename('Lead author' = lead_author,
         'Sample size' = sample_size,
         'Current/former smokers' = current_former_smoker,
         'Never smokers' = never_smoker,
         'Study ID' = study_id)

b$`Lead author` <- to_upper_camel_case(b$`Lead author`, sep_out = ", ")

• Studies reporting partially incomplete smoking data: 13
  These studies typically report smoking history as the measure of exposure which is why they have grouped values for current/former smokers

datatable(b,
          caption = "Studies providing partially complete data on smoking status") %>%
  formatRound(columns = c(2:4), digits = 0, interval = 3, mark = ",")

#Remaining studies
incomplete_smoking_status <- table_1 %>%
  filter(!lead_author %in% full_smoking_status$lead_author) %>%
  filter(!lead_author %in% semi_full_smoking_status$lead_author)

b <- incomplete_smoking_status %>%
  select(lead_author, sample_size, current_smoker, current_former_smoker, former_smoker, never_smoker, never_smoker_unknown, not_stated, missing, study_id) %>%
  rename('Lead author' = lead_author,
         'Sample size' = sample_size,
         'Current smokers' = current_smoker,
         'Current/former smokers' = current_former_smoker,
         'Former smoker' = former_smoker,
         'Never smokers' = never_smoker,
         'Never smoker/unknown' = never_smoker_unknown,
         'Not stated' = not_stated,
         'Missing' = missing,
         'Study ID' = study_id)

b$`Lead author` <- to_upper_camel_case(b$`Lead author`, sep_out = ", ")

• Studies missing key smoking data: 122

datatable(b,
          caption = "Studies providing incomplete data on smoking status") %>%
  formatRound(columns = c(2:8), digits = 0, interval = 3, mark = ",")

0.4 Smoking prevalence by country

#Smoking prevalence by country
country_prevalence_list <- table_1 %>%
    group_by(country, add = T) %>%
  mutate(., current_smok_percentage = current_smoker/total*100) %>%
  mutate(., former_smok_percentage = former_smoker/total*100) %>%
  mutate(., missing_percentage = missing/total*100) %>%
  mutate(., not_stated_percentage = not_stated/total*100) %>%
  mutate(never_smoker_percentage = never_smoker/total*100) %>%
  group_split()

We have split the data into a list where each country can be analysed separately if required

Confidence intervals are obtained through bootstrap estimates of a thousand repetitions of the included studies. Following the production of the confidence intervals the prevalence graphs are produced in the smoking_plots.R script. This uses data extracted on current and former smoking prevalence data from countries which have included studies in this review.

0.4.1 Graph for smoking prevalence

Current smoking prevalence The studies are grouped by the country in which they were conducted. The solid line represents national current smoking prevalence obtained from national smoking survey estimates. Studies are further grouped by the setting in which they were conducted, Hospital (Triangle), Community (Square) and Community and Hospital (Circle). The size of the shape corresponds to the relative size of the study sample. The weighted mean current smoking prevalence is displayed for each setting.

Country smoking prevalence

0.4.2 Graph for former smoking prevalence

Former smoker prevalence This graph is for former smoking prevalence. The data is more sparse for this exposure and therefore several countries are missing for the analysis. Studies are further grouped by the setting in which they were conducted, Hospital (Triangle), Community (Square) and Community and Hospital (Circle). The size of the shape corresponds to the relative size of the study sample. The weighted mean current smoking prevalence is displayed for each setting.

0.5 Final results tables and meta-analyses

0.5.1 Table 1: Included studies

#Updating table 1
table_1_word <- table_1 %>%
  mutate(., current_percentage = current_smoker/total*100,
         former_percentage = former_smoker/total*100,
         current_former_percentage = current_former_smoker/total*100,
         never_smoker_percentage = never_smoker/total*100,
         never_smoker_unknown_percentage = never_smoker_unknown/total*100,
         not_stated_percentage = not_stated/total*100,
         missing_percentage = missing/total*100) %>%
  select(lead_author, date_published, country, sample_size, median_age, iqr_lower, iqr_upper, mean_age, lower_range,
         upper_range, standard_deviation, female_sex_percent, current_percentage, former_percentage, current_former_percentage,
         never_smoker_percentage, never_smoker_unknown_percentage, not_stated_percentage, missing_percentage, study_id) %>%
  replace_na(., list(not_stated_percentage = 0, missing_percentage = 0)) %>%
  mutate(., missing_percentage = not_stated_percentage + missing_percentage,
         missing_percentage = formatC(missing_percentage, digits = 2, format = "f"))

#Smoking completeness
smoking_status <- as_tibble(full_smoking_status$lead_author) %>%
  rename(lead_author = 1) %>%
  mutate(complete_status = c("Yes")) %>%
  bind_rows(., semi_full_smoking_status %>%
              select(lead_author)) %>%
  bind_rows(., incomplete_smoking_status %>%
              select(lead_author)) %>%
  distinct() %>%
  mutate(complete_status = replace_na(complete_status, "No"))

#Data missingness
missingness <- table_6 %>%
  select(lead_author) %>%
  left_join(., table_1_word %>%
              select(lead_author, missing_percentage)) %>%
  mutate(missing_percentage = as.numeric(missing_percentage),
         low_missingness = ifelse(missing_percentage >= 20, "No", "Yes"))

#Defining poor quality
poor_quality <- missingness %>%
  select(-missing_percentage) %>%
  left_join(., smoking_status, by = "lead_author") %>%
  mutate(study_quality = ifelse(low_missingness == "No" | complete_status == "No", "Poor", "Not_poor"))

good_quality <- table_6 %>%
  left_join(., poor_quality, by = "lead_author") %>%
  select(lead_author, study_quality, biochemical_verification, random_representative_sample) %>%
  mutate(study_quality_final = ifelse(study_quality == "Not_poor" & (biochemical_verification == "Yes" | random_representative_sample == "Yes"), "good",
                                      ifelse(study_quality == "Poor", "poor", "fair")))

quality_rating <- good_quality %>%
  select(lead_author, study_quality_final) %>%
  rename("overall_quality" = study_quality_final)

table_1_word <- left_join(table_1_word, good_quality %>%
                            select(lead_author, study_quality_final))


a <- data_study_general %>%
  select(study_id, study_setting)

table_1_word <- left_join(table_1_word, a, by = 'study_id') %>%
  select(1:4, 22, 5:19, 21, 20)
table_1_word$date_published <- as.Date.character(table_1_word$date_published)
write_rds(table_1_word, here::here('data_clean', 'table_1_word.rds'))

table_1_word <- table_1_word %>%
  mutate(median_mean = ifelse(is.na(median_age), mean_age, median_age))


a <- table_1_word %>%
  mutate(median_mean = ifelse(is.na(median_age), mean_age, median_age),
         mean_used = ifelse(is.na(mean_age), '','^'),
         iqr = ifelse(is.na(iqr_lower), NA, paste(iqr_lower, iqr_upper, sep = '-')),
         range_combined = paste(lower_range, upper_range, sep = '-'),
         range_combined = na_if(range_combined, 'NA-NA'),
         standard_deviation = as.numeric(standard_deviation),
         st_dev = paste((as.integer(median_mean-standard_deviation)), as.integer((median_mean+standard_deviation)), sep = '-'),
         st_dev = na_if(st_dev, 'NA-NA'))

a$iqr <- coalesce(a$iqr, a$range_combined, a$st_dev)

a <- a %>%
  select(study_id, lead_author, date_published, country, sample_size, study_setting, median_mean, mean_used, iqr,
         female_sex_percent, current_percentage, current_former_percentage, never_smoker_percentage,
         never_smoker_unknown_percentage, missing_percentage, study_quality_final) %>%
  mutate(median_mean = paste(median_mean, mean_used, sep = ''),
         median_mean = na_if(median_mean, NA)) %>%
  select(-mean_used) %>%
  mutate(median_mean = ifelse(median_mean == 'NA', 'NA', paste(paste(median_mean, iqr, sep = ' ('),')', sep = ''))) %>%
  select(-iqr) %>%
  rename('Study ID' = study_id,
         'Lead author' = lead_author,
         'Date published' = date_published,
         'Country' = country,
         'Sample size' = sample_size,
         'Study setting' = study_setting,
         'Median (IQR)' = median_mean,
         'Female %' = female_sex_percent,
         'Current smoker %' = current_percentage,
         'Current/former smokers %' = current_former_percentage,
         'Never smokers %' = never_smoker_percentage,
         'Never/unknown smokers %' = never_smoker_unknown_percentage,
         'Missing %' = missing_percentage,
         "Study quality" = study_quality_final)
         
a$`Lead author` <- to_upper_camel_case(a$`Lead author`, sep_out = ", ")

cleaned_names <- c('Chow, Us, Cdc',
                   "Gold, Us, Cdc",
                   "Miyara, Updated",
                   "De, La, Rica",
                   "Opensafely, Collaborative",
                   "Bello, Chavolla",
                   "Carrillo, Vega",
                   "De, Lusignan",
                   "Mejia, Vilet",
                   "Heili, Frades",
                   "Vaquero, Roncero",
                   "Al, Hindawi",
                   "Del, Valle",
                   "Soto, Mota",
                   "Eugen, Olsen",
                   "Martinez, Portilla",
                   "Raisi, Estabragh",
                   "Miyara, Medrxiv",
                   "Hernandez, Gardano",
                   "Siso, Almirall",
                   "Martinez, Jiminez",
                   "eugen_olsen",
                   "Mcqueenie",
                   "Martin, Jiminez",
                   "Bello, Chavolla",
                   "De, Melo",
                   "De, Souza",
                   "Bello, Chavolla, Antonio, Villa",
                   "Antonio, Villa")
                   
correct_names <- c('Chow (US CDC)',
                   "Gold (US CDC)",
                   "Miyara",
                   "de la Rica",
                   "The Opensafely Collaborative",
                   "Bello-Chavolla",
                   "Carillo-Vega",
                   "de Lusignan",
                   "Mejia-Vilet",
                   "Heili-Frades",
                   "Vaquero-Roncero",
                   "Al-Hindawi",
                   "del Valle",
                   "Soto-Mota",
                   "Eugen-Olsen",
                   "Martinez-Portilla",
                   "Raisi-Estabragh",
                   "Miyara",
                   "Hernandez-Gardano",
                   "Siso-Almirall",
                   "Martinez-Jiminez",
                   "Eugen-Olsen",
                   "McQueenie",
                   "Martin-Jiminez",
                   "Bello-Chavolla",
                   "de Melo",
                   "de Souza",
                   "Bello-Chavolla, Antonio-Villa",
                   "Antonio-Villa")

author_list <- plyr::mapvalues(a$`Lead author`,
                               from = cleaned_names,
                               to = correct_names)

a$`Lead author` <- author_list

a$`Study setting` <-to_title_case(a$`Study setting`, sep_out = " ")
a$`Country` <-to_title_case(a$`Country`, sep_out = " ")
a$`Country` <- a$Country %>%
  recode('Usa' = 'USA',
         'Uk' = 'UK')

numeric_columns <- c('Median (IQR)', 'Female %', 'Current smoker %', 'Current/former smokers %',
                     'Never smokers %', 'Never/unknown smokers %', 'Missing %')

The total number of included studies is 174. This table is searchable

datatable(a) %>%
  formatRound(
    columns = c(8:14),
    digits = 2,
    interval = 3,
    mark = ",",
    dec.mark = getOption("OutDec")
)

a <- flextable(a) %>%
    set_caption(caption = 'Characteristics of included studies') %>%
  colformat_num(col_keys = numeric_columns, digits = 1, na_str = '-', big.mark = ',') %>%
  colformat_num(col_keys = 'Sample size', digits = 0, na_str = '-', big.mark = ',') %>%
  set_table_properties(width = 1, layout = 'autofit')

save_as_docx(a, path = here('data_clean', 'Table_1.docx'))

0.5.2 Table 2: Smoking and COVID-19 diagnostic testing

#Table 2
table_2_word <-  table_2 %>%
  mutate(., sample = contributing_sample) %>%
  mutate(., negative_test_percentage =formatC(negative_test/sample*100, digits = 2, format = "f"),
         negative_current_percentage = formatC(negative_current_smoker/negative_test*100, digits = 2, format = "f"),
         negative_former_smoker_percentage = formatC(negative_former_smoker/negative_test*100, digits = 2, format = "f"),
         negative_current_former_smoker_percentage = formatC(negative_current_former_smoker/negative_test*100, digits = 2, format = "f"),
         negative_never_smoker_percentage = formatC(negative_never_smoker/negative_test*100, digits = 2, format = "f"),
         negative_not_stated_percentage = formatC(negative_not_stated/negative_test*100, digits = 2, format = "f"),
         positive_test_percentage = formatC(positive_test/sample*100, digits = 2, format = "f"),
         positive_current_smoker_percentage = formatC(positive_current_smoker/positive_test*100, digits = 2, format = "f"),
         positive_former_smoker_percentage = formatC(positive_former_smoker/positive_test*100, digits = 2, format = "f"),
         positive_current_former_smoker_percentage = formatC(positive_current_former_smoker/positive_test*100, digits = 2, format = "f"),
         positive_never_smoker_percentage = formatC(positive_never_smoker/positive_test*100, digits = 2, format = "f"),
         positive_not_stated_percentage = formatC(positive_not_stated/positive_test*100, digits = 2, format = "f")) %>%
  select(-data_on_testing, -missing, -date_published, -sample) %>%
  write_rds(here::here('data_clean', 'table_2_word.rds'))

quality_table_2 <- table_2_word %>%
  left_join(., quality_rating, by = 'lead_author')

a <- table_2_word %>%
  filter(contributing_sample >= 1) %>%
  mutate(Author = lead_author,
         Population_tested = contributing_sample,
         SARS_CoV_2_negative = paste(
           paste(negative_test, 
                 (negative_test_percentage), sep = " ("), "%)", sep = ""),
         N_current_smoker = paste(paste(negative_current_smoker, (negative_current_percentage), sep = " (")
                                  , "%)", sep = ""),
         N_current_smoker = na_if(N_current_smoker, "NA ( NA%)"),
         N_former_smoker = paste(paste(negative_former_smoker, (negative_former_smoker_percentage), sep = " (")
                                 , "%)", sep = ""),
         N_former_smoker = na_if(N_former_smoker, "NA ( NA%)"),
         N_current_former_smoker = paste(paste(negative_current_former_smoker,
                                               (negative_current_former_smoker_percentage), sep = " (")
                                         , "%)", sep = ""),
         N_current_former_smoker = na_if(N_current_former_smoker, "NA ( NA%)"),
         N_never_smoker = paste(paste(negative_never_smoker,
                                      (negative_never_smoker_percentage), sep = " (")
                                , "%)", sep = ""),
         N_never_smoker = na_if(N_never_smoker, "NA ( NA%)"),
         N_not_stated = paste(paste(negative_not_stated,
                                    (negative_not_stated_percentage), sep = " (")
                              , "%)", sep = ""),
         N_not_stated = na_if(N_not_stated, "NA ( NA%)")) %>%
  mutate(SARS_CoV_2_positive = paste(
    paste(positive_test, 
          (positive_test_percentage), sep = " ("), "%)", sep = ""),
    P_current_smoker = paste(paste(positive_current_smoker, (positive_current_smoker_percentage), sep = " (")
                             , "%)", sep = ""),
    P_current_smoker = na_if(P_current_smoker, "NA ( NA%)"),
    P_former_smoker = paste(paste(positive_former_smoker, (positive_former_smoker_percentage), sep = " (")
                            , "%)", sep = ""),
    P_former_smoker = na_if(P_former_smoker, "NA ( NA%)"),
    P_current_former_smoker = paste(paste(positive_current_former_smoker,
                                          (positive_current_former_smoker_percentage), sep = " (")
                                    , "%)", sep = ""),
    P_current_former_smoker = na_if(P_current_former_smoker, "NA ( NA%)"),
    P_never_smoker = paste(paste(positive_never_smoker,
                                 (positive_never_smoker_percentage), sep = " (")
                           , "%)", sep = ""),
    P_never_smoker = na_if(P_never_smoker, "NA ( NA%)"),
    P_not_stated = paste(paste(positive_not_stated,
                               (positive_not_stated_percentage), sep = " (")
                         , "%)", sep = ""),
    P_not_stated = na_if(P_not_stated, "NA ( NA%)")) %>%
  select(Author, Population_tested, SARS_CoV_2_negative, N_current_smoker, N_former_smoker, N_current_former_smoker,
         N_never_smoker, N_not_stated, SARS_CoV_2_positive, P_current_smoker, P_former_smoker, P_current_former_smoker,
         P_never_smoker, P_not_stated)

         
a$Author <- to_upper_camel_case(a$Author, sep_out = ", ")
a$Author  <- a$Author  %>%
  recode("Bello, Chavolla" = "Bello-Chavolla",
         "De, Lusignan" = "de Lusignan",
         "Del, Valle" = "del Valle")

numeric_columns <- c('Population_tested', 'SARS_CoV_2_negative', 'N_current_smoker', 'N_former_smoker',
                     'N_current_former_smoker', 'N_never_smoker', 'N_not_stated', 'SARS_CoV_2_positive', 'P_current_smoker',
                     'P_former_smoker', 'P_current_former_smoker', 'P_never_smoker', 'P_not_stated')
a <- flextable(a) %>%
    set_caption(caption = 'SARS-CoV-2 infection by smoking status') %>%
  colformat_num(j = numeric_columns, digits = 0, na_str = '-', big.mark = ',')

The number of studies reporting on infection by smoking status is: 36

a <- set_header_labels(a,
                       Population_tested = 'Total population tested',
                       SARS_CoV_2_negative = "N (%)",
                       N_current_smoker = "Current smoker (%)",
                       N_former_smoker = "Former smoker (%)",
                       N_current_former_smoker = "Current/former smoker (%)",
                       N_never_smoker = "Never smoker (%)",
                       N_not_stated = "Not stated (%)",
                       SARS_CoV_2_positive = "N (%)",
                       P_current_smoker = "Current smoker (%)",
                       P_former_smoker = "Former smoker (%)",
                       P_current_former_smoker = "Current/former smoker (%)",
                       P_never_smoker = "Never smoker (%)",
                       P_not_stated = "Not stated (%)") %>%
  add_header_row(top = TRUE, values = c("","SARS-CoV-2 negative", "SARS-CoV-2 positive" ), colwidths = c(2, 6, 6)) %>%
  theme_booktabs() %>%
  fix_border_issues() %>%
  set_table_properties(width = 1, layout = 'autofit')  %>%
  colformat_char(na_str = '-')
a

SARS-CoV-2 infection by smoking status

		SARS-CoV-2 negative						SARS-CoV-2 positive
Author	Total population tested	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Not stated (%)	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Not stated (%)
Rentsch	3528	2974 (84.30%)	1444 (48.55%)	704 (23.67%)	-	826 (27.77%)	-	554 (15.70%)	159 (28.70%)	179 (32.31%)	-	216 (38.99%)	-
Fontanet	661	490 (74.13%)	64 (13.06%)	-	-	426 (86.94%)	-	171 (25.87%)	5 (2.92%)	-	-	166 (97.08%)	-
Cho	1331	793 (59.58%)	142 (17.91%)	214 (26.99%)	-	437 (55.11%)	-	538 (40.42%)	111 (20.63%)	145 (26.95%)	-	282 (52.42%)	-
Shah	243	212 (87.24%)	52 (24.53%)	47 (22.17%)	-	113 (53.30%)	-	29 (11.93%)	0 (0.00%)	9 (31.03%)	-	20 (68.97%)	-
Kolin	1474	805 (54.61%)	141 (17.52%)	307 (38.14%)	-	354 (43.98%)	3 (0.37%)	669 (45.39%)	72 (10.76%)	285 (42.60%)	-	303 (45.29%)	9 (1.35%)
de Lusignan	3291	2740 (83.26%)	366 (13.36%)	1450 (52.92%)	-	924 (33.72%)	-	551 (16.74%)	47 (8.53%)	303 (54.99%)	-	201 (36.48%)	-
Valenti	789	689 (87.33%)	197 (28.59%)	-	-	-	492 (71.41%)	40 (5.07%)	7 (17.50%)	-	-	-	33 (82.50%)
Parrotta	76	39 (51.32%)	1 (2.56%)	10 (25.64%)	-	27 (69.23%)	1 (2.56%)	37 (48.68%)	1 (2.70%)	10 (27.03%)	-	25 (67.57%)	1 (2.70%)
Berumen	102875	71353 (69.36%)	-	-	7173 (10.05%)	64180 (89.95%)	-	31522 (30.64%)	-	-	2748 (8.72%)	28774 (91.28%)	-
Israel	24906	20755 (83.33%)	3783 (18.23%)	2671 (12.87%)	-	14301 (68.90%)	-	41151 (165.23%)	406 (0.99%)	483 (1.17%)	-	3262 (7.93%)	-
del Valle	1108	143 (12.91%)	27 (18.88%)	53 (37.06%)	-	-	63 (44.06%)	965 (87.09%)	55 (5.70%)	293 (30.36%)	-	-	617 (63.94%)
Romao	34	20 (58.82%)	-	-	5 (25.00%)	-	15 (75.00%)	14 (41.18%)	-	-	4 (28.57%)	-	10 (71.43%)
Ramlall	11116	4723 (42.49%)	-	-	-	-	-	6393 (57.51%)	-	-	1643.001 (25.70%)	4749.999 (74.30%)	-
Sharma	501	267 (53.29%)	-	-	1 (0.37%)	-	266 (99.63%)	234 (46.71%)	-	-	20 (8.55%)	-	214 (91.45%)
Eugen, Olsen	407	290 (71.25%)	76 (26.21%)	104 (35.86%)	-	102 (35.17%)	-	117 (28.75%)	8 (6.84%)	46 (39.32%)	-	59 (50.43%)	-
Raisi, Estabragh	4510	3184 (70.60%)	-	-	1653 (51.92%)	-	1531 (48.08%)	1326 (29.40%)	-	-	683 (51.51%)	-	643 (48.49%)
Houlihan	177	97 (54.80%)	14 (14.43%)	14 (14.43%)	-	69 (71.13%)	-	80 (45.20%)	7 (8.75%)	19 (23.75%)	-	54 (67.50%)	-
Mcqueenie	428199	424355 (99.10%)	-	-	189299 (44.61%)	235056 (55.39%)	-	1311 (0.31%)	-	-	669 (51.03%)	642 (48.97%)	-
Woolford	4474	3161 (70.65%)	441 (13.95%)	1194 (37.77%)	-	1526 (48.28%)	-	1313 (29.35%)	145 (11.04%)	525 (39.98%)	-	643 (48.97%)	-
Lan	104	83 (79.81%)	-	-	24 (28.92%)	-	59 (71.08%)	21 (20.19%)	-	-	1 (4.76%)	-	20 (95.24%)
Hernandez, Garduno	32583	20279 (62.24%)	-	-	2399 (11.83%)	17861 (88.08%)	-	12304 (37.76%)	-	-	1191 (9.68%)	11083 (90.08%)	-
Govind	6215	6207 (99.87%)	4104 (66.12%)	1669 (26.89%)	-	342 (5.51%)	-	102 (1.64%)	78 (76.47%)	20 (19.61%)	-	2 (1.96%)	-
Gu	4699	3815 (81.19%)	360 (9.44%)	1142 (29.93%)	-	2313 (60.63%)	-	884 (18.81%)	40 (4.52%)	264 (29.86%)	-	580 (65.61%)	-
Kibler	702	680 (96.87%)	25 (3.68%)	-	-	-	655 (96.32%)	22 (3.13%)	1 (4.55%)	-	-	-	21 (95.45%)
Petrilli	10620	5341 (50.29%)	3454 (64.67%)	816 (15.28%)	-	541 (10.13%)	530 (9.92%)	5279 (49.71%)	3268 (61.91%)	902 (17.09%)	-	288 (5.46%)	821 (15.55%)
Bello-Chavolla	150200	98567 (65.62%)	-	-	9624 (9.76%)	-	88943 (90.24%)	51633 (34.38%)	-	-	4366 (8.46%)	-	47267 (91.54%)
Auvinen	61	33 (54.10%)	10 (30.30%)	8 (24.24%)	-	15 (45.45%)	-	28 (45.90%)	1 (3.57%)	9 (32.14%)	-	18 (64.29%)	-
Favara	70	55 (78.57%)	5 (9.09%)	-	-	-	50 (90.91%)	15 (21.43%)	2 (13.33%)	-	-	-	13 (86.67%)
Antonio, Villa	34263	23338 (68.11%)	2293 (9.83%)	-	-	-	21045 (90.17%)	10925 (31.89%)	1023 (9.36%)	-	-	-	9902 (90.64%)
Merzon	7807	7025 (89.98%)	-	-	1136 (16.17%)	-	5889 (83.83%)	782 (10.02%)	-	-	127 (16.24%)	-	655 (83.76%)
Trubiano	2676	2827 (105.64%)	-	-	256 (9.06%)	-	2586 (91.48%)	108 (4.04%)	-	-	3 (2.78%)	-	105 (97.22%)
Shi	1521	1265 (83.17%)	-	-	681 (53.83%)	-	584 (46.17%)	256 (16.83%)	-	-	154 (60.16%)	-	102 (39.84%)
Riley	120620	120461 (99.87%)	2594 (2.15%)	-	-	19914 (16.53%)	97953 (81.32%)	159 (0.13%)	3 (1.89%)	-	-	17 (10.69%)	139 (87.42%)
Alizadehsani	319	196 (61.44%)	-	-	-	-	196 (100.00%)	123 (38.56%)	-	-	1 (0.81%)	-	122 (99.19%)
Merkely	10474	10336 (98.68%)	2904 (28.10%)	2107 (20.39%)	-	5310 (51.37%)	15 (0.15%)	70 (0.67%)	16 (22.86%)	15 (21.43%)	-	38 (54.29%)	1 (1.43%)

save_as_docx(a, path = here('data_clean', 'Table_2.docx'))

0.5.2.1 Smoking and testing meta-analysis

source(here::here('scripts', 'rr_function.R'))

table_2 <- table_2_word

included_studies <- quality_table_2 %>%
  filter(overall_quality != 'poor') %>%
  select(lead_author)

meta <- tibble('author' = table_2$lead_author,
               'negative_smoker' = table_2$negative_current_smoker,
               'negative_never_smoker' = table_2$negative_never_smoker,
               'positive_smoker' = table_2$positive_current_smoker,
               'positive_never_smoker' = table_2$positive_never_smoker,
               'negative_former_smoker' = table_2$negative_former_smoker,
               'positive_former_smoker' = table_2$positive_former_smoker) %>%
        filter(author %in% included_studies$lead_author)
meta$author <- to_upper_camel_case(meta$author, sep_out = ", ")
meta$author <- plyr::mapvalues(meta$author,
                               from = cleaned_names,
                               to = correct_names)

testing_rr <- list()
testing_rr[[1]] <- RR_testing('Rentsch', 'current')
testing_rr[[2]] <- RR_testing('Rentsch', 'former')
testing_rr[[3]] <- RR_testing('Cho', 'current')
testing_rr[[4]] <- RR_testing('Cho', 'former')
testing_rr[[5]] <- RR_testing('Kolin', 'current')
testing_rr[[6]] <- RR_testing('Kolin', 'former')
testing_rr[[7]] <- RR_testing('de Lusignan', 'current')
testing_rr[[8]] <- RR_testing('de Lusignan', 'former')
testing_rr[[9]] <- RR_testing('Parrotta', 'current')
testing_rr[[10]] <- RR_testing('Parrotta', 'former')
testing_rr[[11]] <- RR_testing('Israel', 'current')
testing_rr[[12]] <- RR_testing('Israel', 'former')
testing_rr[[13]] <- RR_testing('Eugen-Olsen', 'current')
testing_rr[[14]] <- RR_testing('Eugen-Olsen', 'former')
testing_rr[[15]] <- RR_testing('Houlihan', 'current')
testing_rr[[16]] <- RR_testing('Houlihan', 'former')
testing_rr[[17]] <- RR_testing('Woolford', 'current')
testing_rr[[18]] <- RR_testing('Woolford', 'former')
testing_rr[[19]] <- RR_testing('Govind', 'current')
testing_rr[[20]] <- RR_testing('Govind', 'former')
testing_rr[[21]] <- RR_testing('Gu', 'current')
testing_rr[[22]] <- RR_testing('Gu', 'former')
testing_rr[[23]] <- RR_testing('Petrilli', 'current')
testing_rr[[24]] <- RR_testing('Petrilli', 'former')
testing_rr[[25]] <- RR_testing('Auvinen', 'current')
testing_rr[[26]] <- RR_testing('Auvinen', 'former')
testing_rr[[27]] <- RR_testing('Merkely', 'current')
testing_rr[[28]] <- RR_testing('Merkely', 'former')

data <- testing_rr
k <- do.call(rbind.data.frame, data)

#SEs for Niedzwiedz et al. 2020

#current vs. never

niedz_log_RR_1<-log(1.15)
niedz_log_SE_1<-(log(1.54)-log(0.86))/3.92

k <- k %>%
  add_row(., study = 'Niedzwiedz', smoking_status = 'current', log_RR = niedz_log_RR_1, log_SE = niedz_log_SE_1)
#former vs. never

niedz_log_RR_2<-log(1.42)
niedz_log_SE_2<-(log(1.69)-log(1.19))/3.92

k <- k %>%
  add_row(., study = 'Niedzwiedz', smoking_status = 'former', log_RR = niedz_log_RR_2, log_SE = niedz_log_SE_2)

numbers_in_analysis <- table_2_word %>%
        left_join(., quality_rating, by = 'lead_author') %>%
        filter(., overall_quality != 'poor') %>%
        add_row(lead_author = 'niedzwiedz', contributing_sample = 1474) %>%
        replace_na(list(negative_not_stated = 0, positive_not_stated = 0)) %>%
        mutate(contributing_sample = (contributing_sample - (negative_not_stated+positive_not_stated))) %>%
                       select(lead_author, contributing_sample) %>%
  na.omit()
numbers_in_analysis$lead_author <- to_upper_camel_case(numbers_in_analysis$lead_author, sep_out = ", ")
numbers_in_analysis$lead_author <- plyr::mapvalues(numbers_in_analysis$lead_author,
                               from = cleaned_names,
                               to = correct_names)

running_meta_count <- k %>%
  group_by(study) %>%
  summarise()

The studies included in meta-analysis are rentsch, cho, kolin, de_lusignan, parrotta, israel, eugen_olsen, houlihan, woolford, govind, gu, petrilli, niedzwiedz, auvinen, and merkely.

• Current versus never smokers

#current vs. never smokers
current_never_meta <- k %>%
  filter(smoking_status == 'current') %>%
  left_join(., numbers_in_analysis %>%
              rename("study" = lead_author),
            by = c("study")) %>%
  mutate(study = paste(study, contributing_sample, sep = ", n = ")) %>%
  select(-contributing_sample)

a <-metagen(current_never_meta$log_RR,
           current_never_meta$log_SE,
           studlab = current_never_meta$study,
           sm="RR",
           comb.fixed = F, comb.random = T)


write_rds(a, here("data_clean", "bayes_testing_current.rds"))

png(here::here('reports', 'figure','fig_3.png'), width=1024, height=546, res=120)
forest(a,
       xlim = c(0.1, 4),
       sortvar = a$TE,
       rightlabs = c('RR', '95% CI', 'Weight'),
       leftlabs = c('Author', 'logRR', 'SE'),
       print.tau2 = F,
       col.diamond = 'blue',
       col.diamond.lines = 'black',
       col.square = 'black',
       col.square.lines = 'black',
       digits.sd = 2,
       colgap.forest.left = unit(15,"mm"))
null <- dev.off()

a

##                           RR            95%-CI %W(random)
## Rentsch, n = 3528     0.4785 [0.3959;  0.5783]        8.4
## Cho, n = 1331         1.1186 [0.9471;  1.3212]        8.5
## Kolin, n = 1462       0.7330 [0.5969;  0.9000]        8.4
## de Lusignan, n = 3291 0.6369 [0.4733;  0.8571]        8.0
## Parrotta, n = 74      1.0400 [0.2528;  4.2784]        2.6
## Israel, n = 24906     0.5218 [0.4734;  0.5753]        8.7
## Eugen-Olsen, n = 407  0.2599 [0.1304;  0.5180]        5.6
## Houlihan, n = 177     0.7593 [0.4016;  1.4356]        6.0
## Woolford, n = 4474    0.8347 [0.7146;  0.9750]        8.5
## Govind, n = 6215      3.2080 [0.7917; 12.9991]        2.7
## Gu, n = 4699          0.4988 [0.3685;  0.6752]        7.9
## Petrilli, n = 9269    1.3994 [1.2707;  1.5412]        8.7
## Auvinen, n = 61       0.1667 [0.0251;  1.1082]        1.7
## Merkely, n = 10458    0.7712 [0.4308;  1.3806]        6.3
## Niedzwiedz, n = 1474  1.1500 [0.8594;  1.5389]        8.0
## 
## Number of studies combined: k = 15
## 
##                          RR           95%-CI     z p-value
## Random effects model 0.7348 [0.5587; 0.9663] -2.21  0.0274
## 
## Quantifying heterogeneity:
##  tau^2 = 0.2222 [0.0863; 0.6632]; tau = 0.4714 [0.2938; 0.8144];
##  I^2 = 95.0% [93.2%; 96.4%]; H = 4.48 [3.83; 5.26]
## 
## Test of heterogeneity:
##       Q d.f.  p-value
##  281.61   14 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Jackson method for confidence interval of tau^2 and tau

Figure 3: Current versus never smokers

• Former versus never smokers

#former vs. never smokers

former_never_meta <- k %>%
  filter(smoking_status == 'former') %>%
  left_join(., numbers_in_analysis %>%
              rename("study" = lead_author),
            by = c("study")) %>%
  mutate(study = paste(study, contributing_sample, sep = ", n = ")) %>%
  select(-contributing_sample)

a<-metagen(former_never_meta$log_RR,
           former_never_meta$log_SE,
           studlab = former_never_meta$study,
           sm="RR", comb.fixed = F, comb.random = T)

write_rds(a, here("data_clean", "bayes_testing_former.rds"))

png(here::here('reports', 'figure','fig_4.png'), width=1024, height=546, res=120)
forest(a,
       xlim = c(0.3, 4),
       sortvar = a$TE,
       rightlabs = c('RR', '95% CI', 'Weight'),
       leftlabs = c('Author', 'logRR', 'SE'),
       print.tau2 = F,
       col.diamond = 'blue',
       col.diamond.lines = 'black',
       col.square = 'black',
       col.square.lines = 'black',
       digits.sd = 2,
       colgap.forest.left = unit(15,"mm"))
null <- dev.off()

a

##                           RR           95%-CI %W(random)
## Rentsch, n = 3528     0.9779 [0.8196; 1.1669]        8.1
## Cho, n = 1331         1.0298 [0.8818; 1.2026]        8.4
## Kolin, n = 1462       1.0439 [0.9281; 1.1741]        9.0
## de Lusignan, n = 3291 0.9674 [0.8229; 1.1373]        8.3
## Parrotta, n = 74      1.0400 [0.6174; 1.7517]        3.5
## Israel, n = 24906     0.8245 [0.7553; 0.9001]        9.3
## Eugen-Olsen, n = 407  0.8368 [0.6108; 1.1465]        5.9
## Houlihan, n = 177     1.3114 [0.9200; 1.8694]        5.3
## Woolford, n = 4474    1.0302 [0.9356; 1.1344]        9.2
## Govind, n = 6215      2.0367 [0.4783; 8.6731]        0.7
## Gu, n = 4699          0.9366 [0.8217; 1.0675]        8.8
## Petrilli, n = 9269    1.5113 [1.3626; 1.6762]        9.1
## Auvinen, n = 61       0.9706 [0.5624; 1.6752]        3.3
## Merkely, n = 10458    0.9948 [0.5484; 1.8046]        2.9
## Niedzwiedz, n = 1474  1.4200 [1.1916; 1.6922]        8.1
## 
## Number of studies combined: k = 15
## 
##                          RR           95%-CI    z p-value
## Random effects model 1.0574 [0.9367; 1.1938] 0.90  0.3668
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0391 [0.0031; 0.0962]; tau = 0.1978 [0.0558; 0.3101];
##  I^2 = 85.6% [77.8%; 90.7%]; H = 2.63 [2.12; 3.27]
## 
## Test of heterogeneity:
##      Q d.f.  p-value
##  97.15   14 < 0.0001
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Jackson method for confidence interval of tau^2 and tau

Figure 4: Former versus never smokers

0.5.3 Table 3: Smoking and COVID-19 hospitalisation

#Table 3
table_3_word <- table_3 %>%
  mutate(., sample = sample_with_outcome) %>%
  mutate(., community_percentage = formatC(number_community/sample*100, digits = 2, format = "f")) %>%
  mutate(., community_current_smoker_percent = formatC(community_current_smoker/number_community*100, digits = 2, format = "f")) %>%
  mutate(., community_former_smoker_percent = formatC(community_former_smoker/number_community*100, digits = 2, format = "f")) %>%
  mutate(., community_current_former_smoker_percent = formatC(community_current_former_smoker/number_community*100, digits = 2, format = "f")) %>%
  mutate(., community_never_smoker_percent = formatC(community_never_smoker/number_community*100, digits = 2, format = "f")) %>%
  mutate(., community_never_unknown_smoker_percent = formatC(community_never_unknown_smoker/number_community*100, digits = 2, format = "f")) %>%
  mutate(., community_not_stated_percent = formatC(community_not_stated/number_community*100, digits = 2, format = "f")) %>%
  mutate(., number_hospitalised_percent = formatC(number_hospitalised/sample*100, digits = 2, format = "f")) %>%
  mutate(., hospitalised_current_smoker_percent = formatC(hospitalised_current_smoker/number_hospitalised*100, digits = 2, format = "f")) %>%
  mutate(., hospitalised_former_smoker_percent = formatC(hospitalised_former_smoker/number_hospitalised*100, digits = 2, format = "f")) %>%
  mutate(., hospitalised_current_former_smoker_percent = formatC(hospitalised_current_former_smoker/number_hospitalised*100, digits = 2, format = "f")) %>%
  mutate(., hospitalised_never_smoker_percent = formatC(hospitalised_never_smoker/number_hospitalised*100, digits = 2, format = "f")) %>%
  mutate(., hospitalised_never_unknown_smoker_percent = formatC(hospitalised_never_unknown_smoker/number_hospitalised*100, digits = 2, format = "f")) %>%
  mutate(., hospitalised_not_stated_percent = formatC(hospitalised_not_stated/number_hospitalised*100, digits = 2, format = "f")) %>%
  select(lead_author, sample_with_outcome, number_community, community_percentage,
         community_current_smoker, community_current_smoker_percent, community_former_smoker, 
         community_former_smoker_percent, community_current_former_smoker,
         community_current_former_smoker_percent, community_never_smoker, community_never_smoker_percent,
         community_never_unknown_smoker, community_never_unknown_smoker_percent, 
         community_not_stated, community_not_stated_percent, number_hospitalised, number_hospitalised_percent,
         hospitalised_current_smoker, hospitalised_current_smoker_percent, hospitalised_former_smoker,
         hospitalised_former_smoker_percent, hospitalised_current_former_smoker, hospitalised_current_former_smoker_percent,
         hospitalised_never_smoker, hospitalised_never_smoker_percent, hospitalised_never_unknown_smoker,
         hospitalised_never_unknown_smoker_percent, hospitalised_not_stated, hospitalised_not_stated_percent) %>%
  write_rds(here::here('data_clean', 'table_3_word.rds'))

The number of studies reporting on hospitalisation by smoking status is: 23

quality_table_3 <- table_3_word %>%
  left_join(., quality_rating, by = 'lead_author') %>%
  select(lead_author, overall_quality)

a <- table_3_word %>%
  filter(number_community >= 1) %>%
  mutate(Author = lead_author,
         Population = sample_with_outcome,
         Community = paste(
           paste(number_community, 
                 as.integer(community_percentage), sep = " ("), "%)", sep = ""),
         C_current_smoker = paste(paste(community_current_smoker, (community_current_smoker_percent), sep = " (")
                                  , "%)", sep = ""),
         C_current_smoker = na_if(C_current_smoker, "NA ( NA%)"),
         C_former_smoker = paste(paste(community_former_smoker, (community_former_smoker_percent), sep = " (")
                                 , "%)", sep = ""),
         C_former_smoker = na_if(C_former_smoker, "NA ( NA%)"),
         C_current_former_smoker = paste(paste(community_current_former_smoker,
                                               (community_current_former_smoker_percent), sep = " (")
                                         , "%)", sep = ""),
         C_current_former_smoker = na_if(C_current_former_smoker, "NA ( NA%)"),
         C_never_smoker = paste(paste(community_never_smoker,
                                      (community_never_smoker_percent), sep = " (")
                                , "%)", sep = ""),
         C_never_smoker = na_if(C_never_smoker, "NA ( NA%)"),
         C_never_unknown_smoker = paste(paste(community_never_unknown_smoker,
                                              (community_never_unknown_smoker_percent), sep = " (")
                                        , "%)", sep = ""),
         C_never_unknown_smoker = na_if(C_never_unknown_smoker, "NA ( NA%)"),
         C_not_stated = paste(paste(community_not_stated,
                                    (community_not_stated_percent), sep = " (")
                              , "%)", sep = ""),
         C_not_stated = na_if(C_not_stated, "NA ( NA%)")) %>%
  mutate(Hospitalised = paste(
    paste(number_hospitalised, 
          as.integer(number_hospitalised_percent), sep = " ("), "%)", sep = ""),
    H_current_smoker = paste(paste(hospitalised_current_smoker, (hospitalised_current_smoker_percent), sep = " (")
                             , "%)", sep = ""),
    H_current_smoker = na_if(H_current_smoker, "NA ( NA%)"),
    H_former_smoker = paste(paste(hospitalised_former_smoker, (hospitalised_former_smoker_percent), sep = " (")
                            , "%)", sep = ""),
    H_former_smoker = na_if(H_former_smoker, "NA ( NA%)"),
    H_current_former_smoker = paste(paste(hospitalised_current_former_smoker,
                                          (hospitalised_current_former_smoker_percent), sep = " (")
                                    , "%)", sep = ""),
    H_current_former_smoker = na_if(H_current_former_smoker, "NA ( NA%)"),
    H_never_smoker = paste(paste(hospitalised_never_smoker,
                                 (hospitalised_never_smoker_percent), sep = " (")
                           , "%)", sep = ""),
    H_never_smoker = na_if(H_never_smoker, "NA ( NA%)"),
    H_never_unknown_smoker = paste(paste(hospitalised_never_unknown_smoker,
                                         (hospitalised_never_unknown_smoker_percent), sep = " (")
                                   , "%)", sep = ""),
    H_never_unknown_smoker = na_if(H_never_unknown_smoker, "NA ( NA%)"),
    H_not_stated = paste(paste(hospitalised_not_stated,
                               (hospitalised_not_stated_percent), sep = " (")
                         , "%)", sep = ""),
    H_not_stated = na_if(H_not_stated, "NA ( NA%)")) %>%
  select(Author, Population, Community, C_current_smoker, C_former_smoker, C_current_former_smoker,
         C_never_smoker,C_never_unknown_smoker,C_not_stated, Hospitalised, H_current_smoker, H_former_smoker,
         H_current_former_smoker, H_never_smoker, H_never_unknown_smoker, H_not_stated)

         
a$Author <- to_upper_camel_case(a$Author, sep_out = ", ")
a$Author <- plyr::mapvalues(a$Author,
                               from = cleaned_names,
                               to = correct_names)

numeric_columns <- 'Population'

a <- flextable(a) %>%
    set_caption(caption = 'COVID-19 hospitalisation by smoking status') %>%
  colformat_num(col_keys = numeric_columns, digits = 0, na_str = '-', big.mark = ',') %>%
  colformat_char(na_str = '-')

a <- set_header_labels(a,
                       Population = 'Population with outcome',
                       Community = "N (%)",
                       C_current_smoker = "Current smoker (%)",
                       C_former_smoker = "Former smoker (%)",
                       C_current_former_smoker = "Current/former smoker (%)",
                       C_never_smoker = "Never smoker (%)",
                       C_never_unknown_smoker = "Never/unknown smoker (%)",
                       C_not_stated = "Not stated (%)",
                       Hospitalised = "N (%)",
                       H_current_smoker = "Current smoker (%)",
                       H_former_smoker = "Former smoker (%)",
                       H_current_former_smoker = "Current/former smoker (%)",
                       H_never_smoker = "Never smoker (%)",
                       H_never_unknown_smoker = "Never/unknown smoker (%)",
                       H_not_stated = "Not stated (%)") %>%
  add_header_row(top = TRUE, values = c("","Community", "Hospitalised" ), colwidths = c(2, 7, 7)) %>%
  theme_booktabs() %>%
  fix_border_issues() %>%
  set_table_properties(width = 1, layout = 'autofit')
a

COVID-19 hospitalisation by smoking status

		Community							Hospitalised
Author	Population with outcome	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Never/unknown smoker (%)	Not stated (%)	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Never/unknown smoker (%)	Not stated (%)
Rentsch	554	269 (48%)	69 (25.65%)	90 (33.46%)	-	110 (40.89%)	-	-	285 (51%)	90 (31.58%)	89 (31.23%)	-	106 (37.19%)	-	-
Chow (US CDC)	6637	5143 (77%)	61 (1.19%)	80 (1.56%)	-	-	-	5002 (97.26%)	1494 (22%)	27 (1.81%)	78 (5.22%)	-	-	-	1389 (92.97%)
Argenziano	1000	151 (15%)	14 (9.27%)	18 (11.92%)	-	119 (78.81%)	-	-	849 (84%)	35 (4.12%)	161 (18.96%)	-	653 (76.91%)	-	-
Lubetzky	54	15 (27%)	-	-	4 (26.67%)	-	-	11 (73.33%)	39 (72%)	-	-	8 (20.51%)	-	-	31 (79.49%)
Carillo-Vega	9946	3922 (39%)	408 (10.40%)	-	-	-	-	3514 (89.60%)	6024 (60%)	486 (8.07%)	-	-	-	-	5538 (91.93%)
Yanover	4353	4180 (96%)	484 (11.58%)	118 (2.82%)	-	3578 (85.60%)	-	-	173 (3%)	30 (17.34%)	11 (6.36%)	-	132 (76.30%)	-	-
Hamer	387109	386349 (99%)	37333 (9.66%)	134542 (34.82%)	-	214474 (55.51%)	-	-	760 (0%)	93 (12.24%)	313 (41.18%)	-	354 (46.58%)	-	-
Heili-Frades	4712	1973 (41%)	121 (6.13%)	222 (11.25%)	-	-	1630 (82.62%)	1630 (82.62%)	2739 (58%)	112 (4.09%)	598 (21.83%)	-	-	2029 (74.08%)	-
Freites	123	69 (56%)	1 (1.45%)	-	-	-	-	68 (98.55%)	54 (43%)	3 (5.56%)	-	-	-	-	51 (94.44%)
Berumen	102875	18832 (18%)	-	-	1546 (8.21%)	-	17286 (91.79%)	-	12690 (12%)	-	-	1202 (9.47%)	-	11488 (90.53%)	-
Gianfrancesco	600	323 (53%)	-	-	61 (18.89%)	-	-	262 (81.11%)	277 (46%)	-	-	68 (24.55%)	-	-	209 (75.45%)
Chaudhry	40	19 (47%)	-	-	0 (0.00%)	-	-	19 (100.00%)	21 (52%)	-	-	6 (28.57%)	-	-	15 (71.43%)
Giannouchos	89756	58485 (65%)	4679 (8.00%)	-	-	-	53806 (92.00%)	-	31271 (34%)	2721 (8.70%)	-	-	-	28550 (91.30%)	-
Wang, Oekelen	57	22 (38%)	-	-	6 (27.27%)	-	-	16 (72.73%)	36 (63%)	-	-	15 (41.67%)	-	-	20 (55.56%)
Miyara	470	132 (28%)	14 (10.61%)	41 (31.06%)	-	77 (58.33%)	-	-	338 (71%)	18 (5.33%)	111 (32.84%)	-	209 (61.83%)	-	-
Suleyman	463	108 (23%)	-	-	23 (21.30%)	-	-	85 (78.70%)	355 (76%)	-	-	137 (38.59%)	-	-	218 (61.41%)
Garassino	196	48 (24%)	10 (20.83%)	27 (56.25%)	-	11 (22.92%)	-	-	152 (77%)	38 (25.00%)	84 (55.26%)	-	26 (17.11%)	-	-
Siso-Almirall	260	119 (45%)	-	-	31 (26.05%)	-	-	88 (73.95%)	141 (54%)	-	-	50 (35.46%)	-	-	91 (64.54%)
Gu	884	511 (57%)	30 (5.87%)	126 (24.66%)	-	355 (69.47%)	-	-	373 (42%)	10 (2.68%)	138 (37.00%)	-	225 (60.32%)	-	-
Killerby	531	311 (58%)	-	-	37 (11.90%)	222 (71.38%)	-	52 (16.72%)	220 (41%)	-	-	54 (24.55%)	157 (71.36%)	-	9 (4.09%)
Petrilli	5279	2538 (48%)	147 (5.79%)	337 (13.28%)	-	1678 (66.12%)	-	376 (14.81%)	2741 (51%)	141 (5.14%)	565 (20.61%)	-	1590 (58.01%)	-	445 (16.23%)
Nguyen	689	333 (48%)	-	-	57 (17.12%)	-	-	276 (82.88%)	356 (51%)	-	-	114 (32.02%)	-	-	242 (67.98%)
Mendy	689	473 (68%)	-	-	84 (17.76%)	-	-	389 (82.24%)	216 (31%)	-	-	86 (39.81%)	-	-	130 (60.19%)

save_as_docx(a, path = here('data_clean', 'Table_3.docx'))

0.5.3.1 Smoking and hospitalisation meta-analysis

The studies included in meta-analysis are rentsch, argenziano, yanover, hamer, miyara_medrxiv, garassino, gu, and petrilli.

# Data --------------------------------------------------------------------
meta <- tibble('author' = table_3$lead_author,
               'community_smoker' = table_3$community_current_smoker,
               'community_never_smoker' = table_3$community_never_smoker, 
               'hospitalised_smoker' = table_3$hospitalised_current_smoker, 
               'hospitalised_never_smoker' = table_3$hospitalised_never_smoker,
               'community_former_smoker' = table_3$community_former_smoker,
               'hospitalised_former_smoker' = table_3$hospitalised_former_smoker) %>%
  filter(author %in% included_studies$lead_author)

meta$author <- to_upper_camel_case(meta$author, sep_out = ", ")
meta$author <- plyr::mapvalues(meta$author,
                               from = cleaned_names,
                               to = correct_names)  

# Current smoker hospitalisation ------------------------------------------
event_rates_smoker <- meta %>%
  mutate(., Ee = hospitalised_smoker) %>%
  mutate(., Ne = (hospitalised_smoker+community_smoker)) %>%
  mutate(., Ec = hospitalised_never_smoker) %>%
  mutate(., Nc = (hospitalised_never_smoker+community_never_smoker)) %>%
  rename('Author' = author) %>%
  select(Author, Ee, Ne, Ec, Nc)

event_rates_smoker <- metabin(Ee,
                              Ne,
                              Ec,
                              Nc,
                              data = event_rates_smoker,
                              studlab = paste(Author),
                              comb.fixed = F,
                              comb.random = T,
                              method.tau = 'SJ',
                              hakn = T,
                              prediction = F,
                              incr = 0.1,
                              sm = 'RR')

write_rds(event_rates_smoker, here("data_clean", "bayes_hospital_current.rds"))

png(here::here('reports', 'figure', 'fig_5.png'), width=1480, height=546, res=120)
current_smoker_hospitalisation <- forest(event_rates_smoker,
                                sortvar = Author,
                                xlim = c(0.3, 3),
                                rightlabs = c('RR', '95% CI', 'Weight'),
                                leftlabs = c('Author', 'logRR', 'SE'),
                                lab.e = 'Current smoker',
                                lab.c = 'Never smoker',
                                print.tau2 = F,
                                col.diamond = 'blue',
                                col.diamond.lines = 'black',
                                col.square = 'black',
                                col.square.lines = 'black',
                                digits.sd = 2)
null <- dev.off()

event_rates_smoker

##                RR           95%-CI %W(random)
## Rentsch    1.1534 [0.9517; 1.3980]       14.1
## Argenziano 0.8445 [0.7056; 1.0106]       14.3
## Yanover    1.6404 [1.1156; 2.4121]       10.4
## Hamer      1.5080 [1.2004; 1.8944]       13.5
## Miyara     0.7697 [0.5626; 1.0532]       11.8
## Garassino  1.1266 [0.8731; 1.4537]       13.0
## Gu         0.6444 [0.3732; 1.1130]        7.8
## Petrilli   1.0063 [0.8897; 1.1380]       15.1
## 
## Number of studies combined: k = 8
## 
##                          RR           95%-CI    t p-value
## Random effects model 1.0558 [0.8236; 1.3535] 0.52  0.6209
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0735 [0.0146; 0.3895]; tau = 0.2711 [0.1210; 0.6241];
##  I^2 = 75.8% [51.6%; 87.9%]; H = 2.03 [1.44; 2.88]
## 
## Test of heterogeneity:
##      Q d.f. p-value
##  28.96    7  0.0001
## 
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - Sidik-Jonkman estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
## - Hartung-Knapp adjustment for random effects model

Figure 5: Current versus never smokers

# Former smoker hospitalisation -------------------------------------------
event_rates_former <- meta %>%
  mutate(., Ee = hospitalised_former_smoker) %>%
  mutate(., Ne = (hospitalised_former_smoker+community_former_smoker)) %>%
  mutate(., Ec = hospitalised_never_smoker) %>%
  mutate(., Nc = (hospitalised_never_smoker+community_never_smoker)) %>%
  rename('Author' = author) %>%
  select(Author, Ee, Ne, Ec, Nc)

event_rates_former <- metabin(Ee,
                              Ne,
                              Ec,
                              Nc,
                              data = event_rates_former,
                              studlab = paste(Author),
                              comb.fixed = F,
                              comb.random = T,
                              method.tau = 'SJ',
                              hakn = T,
                              prediction = F,
                              incr = 0.1,
                              sm = 'RR')

write_rds(event_rates_former, here("data_clean", "bayes_hospital_former.rds"))

png(here::here('reports', 'figure', 'fig_6.png'), width=1480, height=546, res=120)
former_smoker_hospitalisation <- forest(event_rates_former,
                                         sortvar = Author,
                                         xlim = c(0.3, 3),
                                         rightlabs = c('RR', '95% CI', 'Weight'),
                                         leftlabs = c('Author', 'logRR', 'SE'),
                                         lab.e = 'Former smoker',
                                         lab.c = 'Never smoker',
                                         print.tau2 = F,
                                         col.diamond = 'blue',
                                         col.diamond.lines = 'black',
                                         col.square = 'black',
                                         col.square.lines = 'black',
                                         digits.sd = 2)
null <- dev.off()
event_rates_former

##                RR           95%-CI %W(random)
## Rentsch    1.0132 [0.8292; 1.2380]       12.4
## Argenziano 1.0634 [1.0039; 1.1263]       14.9
## Yanover    2.3966 [1.3292; 4.3214]        5.2
## Hamer      1.4085 [1.2100; 1.6396]       13.5
## Miyara     0.9993 [0.8867; 1.1262]       14.1
## Garassino  1.0769 [0.8517; 1.3617]       11.7
## Gu         1.3475 [1.1551; 1.5719]       13.4
## Petrilli   1.2874 [1.2107; 1.3691]       14.8
## 
## Number of studies combined: k = 8
## 
##                          RR           95%-CI    t p-value
## Random effects model 1.2077 [1.0038; 1.4531] 2.41  0.0466
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0478 [0.0065; 0.2847]; tau = 0.2186 [0.0808; 0.5336];
##  I^2 = 84.0% [70.0%; 91.4%]; H = 2.50 [1.83; 3.41]
## 
## Test of heterogeneity:
##      Q d.f.  p-value
##  43.63    7 < 0.0001
## 
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - Sidik-Jonkman estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
## - Hartung-Knapp adjustment for random effects model

running_meta_count <- running_meta_count %>%
  full_join(., meta %>%
  rename("study" = author), by = "study") %>%
  select(study)

Figure 6: Former versus never smokers

0.5.4 Table 4: Smoking and COVID-19 disease severity

#Table 4
table_4_word <- table_4 %>%
  mutate(., sample = sample_with_severity) %>%
  mutate(., non_severe_disease_percentage = formatC(non_severe_disease/sample*100, digits = 2, format = "f")) %>%
  mutate(., non_severe_current_smoker_percent = formatC(non_severe_current_smoker/non_severe_disease*100, digits = 2, format = "f")) %>%
  mutate(., non_severe_former_smoker_percent = formatC(non_severe_former_smoker/non_severe_disease*100, digits = 2, format = "f")) %>%
  mutate(., non_severe_current_former_smoker_percent = formatC(non_severe_current_former_smoker/non_severe_disease*100, digits = 2, format = "f")) %>%
  mutate(., non_severe_never_smoker_percent = formatC(non_severe_never_smoker/non_severe_disease*100, digits = 2, format = "f")) %>%
  mutate(., non_severe_never_unknown_smoker_percent = formatC(non_severe_never_unknown_smoker/non_severe_disease*100, digits = 2, format = "f")) %>%
  mutate(., non_severe_not_stated_percent = formatC(non_severe_not_stated/non_severe_disease*100, digits = 2, format = "f")) %>%
  mutate(., severe_disease_number_percent = formatC(severe_disease_number/sample*100, digits = 2, format = "f")) %>%
  mutate(., severe_disease_current_smoker_percent = formatC(severe_disease_current_smoker/severe_disease_number*100, digits = 2, format = "f")) %>%
  mutate(., severe_disease_former_smoker_percent = formatC(severe_disease_former_smoker/severe_disease_number*100, digits = 2, format = "f")) %>%
  mutate(., severe_disease_current_former_smoker_percent = formatC(severe_disease_current_former_smoker/severe_disease_number*100, digits = 2, format = "f")) %>%
  mutate(., severe_disease_never_smoker_percent = formatC(severe_disease_never_smoker/severe_disease_number*100, digits = 2, format = "f")) %>%
  mutate(., severe_disease_never_unknown_percent = formatC(severe_disease_never_unknown/severe_disease_number*100, digits = 2, format = "f")) %>%
  mutate(., severe_disease_not_stated_percent = formatC(severe_disease_not_stated/severe_disease_number*100, digits = 2, format = "f")) %>%
  select(lead_author, sample, non_severe_disease, non_severe_disease_percentage, non_severe_current_smoker,
         non_severe_current_smoker_percent, non_severe_former_smoker, non_severe_former_smoker_percent,
         non_severe_current_former_smoker, non_severe_current_former_smoker_percent, non_severe_never_smoker,
         non_severe_never_smoker_percent, non_severe_never_unknown_smoker, non_severe_never_unknown_smoker_percent,
         non_severe_not_stated, non_severe_not_stated_percent, severe_disease_number, severe_disease_number_percent,
         severe_disease_current_smoker, severe_disease_current_smoker_percent, 
         severe_disease_former_smoker, severe_disease_former_smoker_percent, severe_disease_current_former_smoker,
         severe_disease_current_former_smoker_percent, severe_disease_never_smoker, severe_disease_never_smoker_percent,
         severe_disease_never_unknown, severe_disease_never_unknown_percent, severe_disease_not_stated,
         severe_disease_not_stated_percent) %>%
  write_rds(here::here('data_clean', 'table_4_word.rds'))

The number of studies reporting on hospitalisation by smoking status is: 44

quality_table_4 <- table_4_word %>%
  left_join(., quality_rating, by = 'lead_author') %>%
  select(lead_author, overall_quality)

a <- table_4_word %>%
  mutate(Author = lead_author,
         Population = sample,
         non_severe_disease = paste(
           paste(non_severe_disease, 
                 as.integer(non_severe_disease_percentage), sep = " ("), "%)", sep = ""),
         ns_current_smoker = paste(paste(non_severe_current_smoker, (non_severe_current_smoker_percent), sep = " (")
                                  , "%)", sep = ""),
         ns_current_smoker = na_if(ns_current_smoker, "NA ( NA%)"),
         ns_former_smoker = paste(paste(non_severe_former_smoker, (non_severe_former_smoker_percent), sep = " (")
                                 , "%)", sep = ""),
         ns_former_smoker = na_if(ns_former_smoker, "NA ( NA%)"),
         ns_current_former_smoker = paste(paste(non_severe_current_former_smoker,
                                               (non_severe_current_former_smoker_percent), sep = " (")
                                         , "%)", sep = ""),
         ns_current_former_smoker = na_if(ns_current_former_smoker, "NA ( NA%)"),
         ns_never_smoker = paste(paste(non_severe_never_smoker,
                                      (non_severe_never_smoker_percent), sep = " (")
                                , "%)", sep = ""),
         ns_never_smoker = na_if(ns_never_smoker, "NA ( NA%)"),
         ns_never_unknown_smoker = paste(paste(non_severe_never_unknown_smoker,
                                              (non_severe_never_unknown_smoker_percent), sep = " (")
                                        , "%)", sep = ""),
         ns_never_unknown_smoker = na_if(ns_never_unknown_smoker, "NA ( NA%)"),
         ns_not_stated = paste(paste(non_severe_not_stated,
                                    (non_severe_not_stated_percent), sep = " (")
                              , "%)", sep = ""),
         ns_not_stated = na_if(ns_not_stated, "NA ( NA%)"),
         ) %>%
  mutate(severe_disease = paste(
    paste(severe_disease_number, 
          as.integer(severe_disease_number_percent), sep = " ("), "%)", sep = ""),
    s_current_smoker = paste(paste(severe_disease_current_smoker, (severe_disease_current_smoker_percent), sep = " (")
                             , "%)", sep = ""),
    s_current_smoker = na_if(s_current_smoker, "NA ( NA%)"),
    s_former_smoker = paste(paste(severe_disease_former_smoker, (severe_disease_former_smoker_percent), sep = " (")
                            , "%)", sep = ""),
    s_former_smoker = na_if(s_former_smoker, "NA ( NA%)"),
    s_current_former_smoker = paste(paste(severe_disease_current_former_smoker,
                                          (severe_disease_current_former_smoker_percent), sep = " (")
                                    , "%)", sep = ""),
    s_current_former_smoker = na_if(s_current_former_smoker, "NA ( NA%)"),
    s_never_smoker = paste(paste(severe_disease_never_smoker,
                                 (severe_disease_never_smoker_percent), sep = " (")
                           , "%)", sep = ""),
    s_never_smoker = na_if(s_never_smoker, "NA ( NA%)"),
    s_never_unknown_smoker = paste(paste(severe_disease_never_unknown,
                                         (severe_disease_never_unknown_percent), sep = " (")
                                   , "%)", sep = ""),
    s_never_unknown_smoker = na_if(s_never_unknown_smoker, "NA ( NA%)"),
    s_not_stated = paste(paste(severe_disease_not_stated,
                                 (severe_disease_not_stated_percent), sep = " (")
                           , "%)", sep = ""),
    s_not_stated = na_if(s_not_stated, "NA ( NA%)")) %>%
  select(Author, Population, non_severe_disease, ns_current_smoker, ns_former_smoker,
         ns_current_former_smoker, ns_never_smoker, ns_never_unknown_smoker, ns_not_stated,
         severe_disease, s_current_smoker, s_former_smoker, s_current_former_smoker, s_never_smoker, s_never_unknown_smoker,
         s_not_stated)

         
a$Author <- to_upper_camel_case(a$Author, sep_out = ", ")
a$Author <- plyr::mapvalues(a$Author,
                               from = cleaned_names,
                               to = correct_names)  
numeric_columns <- 'Population'

a <- flextable(a) %>%
    set_caption(caption = 'COVID-19 severity by smoking status') %>%
  colformat_num(col_keys = numeric_columns, digits = 0, na_str = '-', big.mark = ',') %>%
  colformat_char(na_str = '-')

a <- set_header_labels(a,
                       Population = "Population with severity",
                       non_severe_disease = "N (%)",
                       ns_current_smoker = "Current smoker (%)",
                       ns_former_smoker = "Former smoker (%)",
                       ns_current_former_smoker = "Current/former smoker (%)",
                       ns_never_unknown_smoker = "Never/unknown smoker (%)",
                       ns_never_smoker = "Never smoker (%)",
                       ns_not_stated = "Not stated (%)",
                       severe_disease = "N (%)", 
                       s_current_smoker = "Current smoker (%)", 
                       s_former_smoker = "Former smoker (%)", 
                       s_current_former_smoker = "Current/former smoker (%)",
                       s_never_unknown_smoker = "Never/unknown smoker (%)",
                       s_never_smoker = "Never smoker (%)",
                       s_not_stated = "Not stated (%)") %>%
  add_header_row(top = TRUE, values = c("","Non severe disease", "Severe disease" ), colwidths = c(2, 7, 7)) %>%
  theme_booktabs() %>%
  fix_border_issues() %>%
  set_table_properties(width = 1, layout = 'autofit')
a

COVID-19 severity by smoking status

		Non severe disease							Severe disease
Author	Population with severity	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Never/unknown smoker (%)	Not stated (%)	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Never/unknown smoker (%)	Not stated (%)
Guan, Ni	1085	913 (84%)	108 (11.83%)	12 (1.31%)	-	793 (86.86%)	-	-	172 (15%)	29 (16.86%)	9 (5.23%)	-	134 (77.91%)	-	-
Zhang, Dong	9	3 (33%)	0 (0.00%)	3 (100.00%)	-	0 (0.00%)	-	-	6 (66%)	2 (33.33%)	4 (66.67%)	-	0 (0.00%)	-	-
Wan	9	8 (88%)	8 (100.00%)	0 (0.00%)	-	0 (0.00%)	-	-	1 (11%)	1 (100.00%)	0 (0.00%)	-	0 (0.00%)	-	-
Huang, Wang	3	3 (100%)	3 (100.00%)	0 (0.00%)	-	0 (0.00%)	-	-	0 (0%)	0 (NaN%)	0 (NaN%)	-	0 (NaN%)	-	-
Rentsch	285	168 (58%)	47 (27.98%)	53 (31.55%)	-	68 (40.48%)	-	-	117 (41%)	43 (36.75%)	36 (30.77%)	-	38 (32.48%)	-	-
Hu	323	151 (46%)	-	-	12 (7.95%)	-	139 (92.05%)	-	172 (53%)	-	-	26 (15.12%)	-	146 (84.88%)	-
Wang, Pan	125	100 (80%)	-	-	9 (9.00%)	-	91 (91.00%)	-	25 (20%)	-	-	7 (28.00%)	-	18 (72.00%)	-
Kim	27	21 (77%)	3 (14.29%)	-	-	-	18 (85.71%)	-	6 (22%)	2 (33.33%)	0 (0.00%)	-	-	4 (66.67%)	-
Shi, Yu	474	425 (89%)	-	-	34 (8.00%)	-	391 (92.00%)	-	49 (10%)	-	-	6 (12.24%)	-	43 (87.76%)	-
Liao, Feng	148	92 (62%)	-	-	5 (5.43%)	-	-	87 (94.57%)	56 (37%)	3 (5.36%)	-	-	-	-	53 (94.64%)
Shi, Ren	134	88 (65%)	-	-	8 (9.09%)	-	-	80 (90.91%)	46 (34%)	-	-	6 (13.04%)	-	-	40 (86.96%)
Hadjadj	50	15 (30%)	1 (6.67%)	2 (13.33%)	-	12 (80.00%)	-	-	35 (70%)	0 (0.00%)	7 (20.00%)	-	28 (80.00%)	-	-
Zheng, Xiong	73	43 (58%)	-	-	6 (13.95%)	37 (86.05%)	-	-	30 (41%)	-	-	2 (6.67%)	28 (93.33%)	-	-
de la Rica	48	26 (54%)	-	-	6 (23.08%)	-	-	20 (76.92%)	20 (41%)	-	-	4 (20.00%)	-	-	16 (80.00%)
Yin, Yang	106	47 (44%)	-	-	6 (12.77%)	-	-	41 (87.23%)	59 (55%)	-	-	12 (20.34%)	-	-	47 (79.66%)
Allenbach	147	100 (68%)	-	-	9 (9.00%)	-	-	91 (91.00%)	47 (31%)	-	-	0 (0.00%)	-	-	47 (100.00%)
Goyal	393	263 (66%)	14 (5.32%)	-	-	-	-	249 (94.68%)	130 (33%)	6 (4.62%)	-	-	-	-	124 (95.38%)
Feng	454	333 (73%)	27 (8.11%)	-	-	-	-	306 (91.89%)	121 (26%)	17 (14.05%)	-	-	-	-	104 (85.95%)
Yao	108	83 (76%)	1 (1.20%)	-	-	-	-	82 (98.80%)	25 (23%)	3 (12.00%)	-	-	-	-	22 (88.00%)
Sami	490	400 (81%)	53 (13.25%)	-	-	-	-	347 (86.75%)	90 (18%)	16 (17.78%)	-	-	-	-	74 (82.22%)
Regina	200	163 (81%)	9 (5.52%)	-	-	-	-	154 (94.48%)	37 (18%)	0 (0.00%)	-	-	-	-	37 (100.00%)
Feuth	28	21 (75%)	1 (4.76%)	7 (33.33%)	-	13 (61.90%)	-	-	7 (25%)	2 (28.57%)	1 (14.29%)	-	4 (57.14%)	-	-
Mejia-Vilet	329	214 (65%)	-	-	13 (6.07%)	-	-	201 (93.93%)	115 (34%)	-	-	10 (8.70%)	-	-	105 (91.30%)
Chen, Jiang	135	54 (40%)	-	-	4 (7.41%)	-	-	50 (92.59%)	81 (60%)	-	-	9 (11.11%)	-	-	72 (88.89%)
Vaquero-Roncero	146	75 (51%)	-	-	4 (5.33%)	-	-	71 (94.67%)	71 (48%)	-	-	6 (8.45%)	-	-	65 (91.55%)
Kim, Garg	2490	1692 (67%)	112 (6.62%)	395 (23.35%)	-	-	1185 (70.04%)	-	798 (32%)	38 (4.76%)	247 (30.95%)	-	-	512 (64.16%)	-
Wu	174	92 (52%)	-	-	47 (51.09%)	-	45 (48.91%)	-	82 (47%)	11 (13.41%)	-	-	-	71 (86.59%)	-
Chaudhry	40	34 (85%)	-	-	5 (14.71%)	-	-	29 (85.29%)	6 (15%)	-	-	1 (16.67%)	-	-	5 (83.33%)
Garibaldi	832	532 (63%)	25 (4.70%)	107 (20.11%)	-	-	-	400 (75.19%)	300 (36%)	21 (7.00%)	81 (27.00%)	-	-	-	198 (66.00%)
Kuderer	928	686 (73%)	35 (5.10%)	210 (30.61%)	-	370 (53.94%)	-	29 (4.23%)	242 (26%)	8 (3.31%)	116 (47.93%)	-	99 (40.91%)	15 (6.20%)	4 (1.65%)
Romao	14	14 (100%)	-	-	4 (28.57%)	-	-	10 (71.43%)	0 (0%)	-	-	-	-	-	-
Giannouchos	89756	78050 (86%)	6322 (8.10%)	-	-	-	71728 (91.90%)	-	11706 (13%)	1089 (9.30%)	-	-	-	10617 (90.70%)	-
Cen	1007	720 (71%)	-	-	70 (9.72%)	-	-	650 (90.28%)	287 (28%)	-	-	18 (6.27%)	-	-	269 (93.73%)
Maraschini	132	89 (67%)	-	11 (12.36%)	-	78 (87.64%)	-	-	43 (32%)	-	3 (6.98%)	-	40 (93.02%)	-	-
Siso-Almirall	260	212 (81%)	-	-	60 (28.30%)	-	-	152 (71.70%)	48 (18%)	-	-	21 (43.75%)	-	-	27 (56.25%)
Gu	884	511 (57%)	30 (5.87%)	126 (24.66%)	-	355 (69.47%)	-	-	134 (15%)	3 (2.24%)	61 (45.52%)	-	70 (52.24%)	-	-
Petrilli	2729	1739 (63%)	97 (5.58%)	325 (18.69%)	-	1067 (61.36%)	-	250 (14.38%)	990 (36%)	44 (4.44%)	236 (23.84%)	-	517 (52.22%)	-	193 (19.49%)
Mendy	689	598 (86%)	-	-	133 (22.24%)	-	-	465 (77.76%)	91 (13%)	-	-	37 (40.66%)	-	-	54 (59.34%)
Pongpirul	193	161 (83%)	-	-	25 (15.53%)	106 (65.84%)	-	30 (18.63%)	32 (16%)	-	-	4 (12.50%)	21 (65.62%)	-	7 (21.88%)
Jin, Gu	6	2 (33%)	-	-	0 (0.00%)	-	-	4 (200.00%)	4 (66%)	-	-	2 (50.00%)	-	-	2 (50.00%)
Senkal	611	446 (73%)	48 (10.76%)	-	-	-	-	398 (89.24%)	165 (27%)	21 (12.73%)	-	-	-	-	144 (87.27%)
Patel	129	89 (68%)	26 (29.21%)	-	-	-	58 (65.17%)	5 (5.62%)	40 (31%)	22 (55.00%)	-	-	-	14 (35.00%)	4 (10.00%)
Maucourant	27	10 (37%)	1 (10.00%)	2 (20.00%)	-	2 (20.00%)	-	5 (50.00%)	17 (62%)	2 (11.76%)	5 (29.41%)	-	9 (52.94%)	-	1 (5.88%)
Xie	619	469 (75%)	-	-	32 (6.82%)	-	-	437 (93.18%)	150 (24%)	-	-	19 (12.67%)	-	-	131 (87.33%)

save_as_docx(a, path = here('data_clean', 'Table_4.docx'))

0.5.4.1 Smoking and severity meta-analysis

The studies included in meta-analysis are guan_ni, rentsch, hadjadj, feuth, kuderer, gu, and petrilli.

# Data --------------------------------------------------------------------
meta <- tibble('author' = table_4$lead_author,
               'non_severe_smoker' = table_4$non_severe_current_smoker,
               'non_severe_never_smoker' = table_4$non_severe_never_smoker,
               'severe_smoker' = table_4$severe_disease_current_smoker,
               'severe_never_smoker' = table_4$severe_disease_never_smoker,
               'non_severe_former_smoker' = table_4$non_severe_former_smoker,
               'severe_former_smoker' = table_4$severe_disease_former_smoker) %>%
  filter(author %in% included_studies$lead_author)

meta$author <- to_upper_camel_case(meta$author, sep_out = ", ")
meta$author <- plyr::mapvalues(meta$author,
                               from = cleaned_names,
                               to = correct_names)  

# Current smoker severity ------------------------------------------
event_rates_smoker <- meta %>%
  mutate(., Ee = severe_smoker) %>%
  mutate(., Ne = (severe_smoker+non_severe_smoker)) %>%
  mutate(., Ec = severe_never_smoker) %>%
  mutate(., Nc = (severe_never_smoker+non_severe_never_smoker)) %>%
  rename('Author' = author) %>%
  select(Author, Ee, Ne, Ec, Nc)

event_rates_smoker <- metabin(Ee,
                              Ne,
                              Ec,
                              Nc,
                              data = event_rates_smoker,
                              studlab = paste(Author),
                              comb.fixed = F,
                              comb.random = T,
                              method.tau = 'SJ',
                              hakn = F,
                              prediction = F,
                              incr = 0.1,
                              sm = 'RR')

write_rds(event_rates_smoker, here("data_clean", "bayes_severity_current.rds"))

png(here::here('reports', 'figure', 'fig_7.png'), width=1024, height=546, res=120)
current_smoker_severity <- forest(event_rates_smoker,
                                         sortvar = Author,
                                         xlim = c(0.1, 5),
                                         rightlabs = c('RR', '95% CI', 'Weight'),
                                         leftlabs = c('Author', 'logRR', 'SE'),
                                         lab.e = 'Current smoker',
                                         lab.c = 'Never smoker',
                                         print.tau2 = F,
                                         col.diamond = 'blue',
                                         col.diamond.lines = 'black',
                                         col.square = 'black',
                                         col.square.lines = 'black',
                                         digits.sd = 2)
null <- dev.off()

event_rates_smoker

##              RR            95%-CI %W(random)
## Guan, Ni 1.4644 [1.0226;  2.0970]       20.5
## Rentsch  1.3327 [0.9544;  1.8611]       20.9
## Hadjadj  0.1297 [0.0004; 47.9753]        0.6
## Feuth    2.8333 [0.8772;  9.1514]        9.6
## Kuderer  0.8814 [0.4605;  1.6870]       16.2
## Gu       0.5519 [0.1837;  1.6581]       10.3
## Petrilli 0.9561 [0.7408;  1.2339]       21.8
## 
## Number of studies combined: k = 7
## 
##                          RR           95%-CI    z p-value
## Random effects model 1.1432 [0.7112; 1.8377] 0.55  0.5804
## 
## Quantifying heterogeneity:
##  tau^2 = 0.2521 [0.0000; 1.6188]; tau = 0.5021 [0.0000; 1.2723];
##  I^2 = 37.8% [0.0%; 73.8%]; H = 1.27 [1.00; 1.95]
## 
## Test of heterogeneity:
##     Q d.f. p-value
##  9.64    6  0.1404
## 
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - Sidik-Jonkman estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
## - Continuity correction of 0.1 in studies with zero cell frequencies

Figure 7: Current versus never smokers

# Former smoker hospitalisation -------------------------------------------
event_rates_former <- meta %>%
  mutate(., Ee = severe_former_smoker) %>%
  mutate(., Ne = (severe_former_smoker+non_severe_former_smoker)) %>%
  mutate(., Ec = severe_never_smoker) %>%
  mutate(., Nc = (severe_never_smoker+non_severe_never_smoker)) %>%
  rename('Author' = author) %>%
  select(Author, Ee, Ne, Ec, Nc)

event_rates_former <- metabin(Ee,
                              Ne,
                              Ec,
                              Nc,
                              data = event_rates_former,
                              studlab = paste(Author),
                              comb.fixed = F,
                              comb.random = T,
                              method.tau = 'SJ',
                              hakn = F,
                              prediction = F,
                              incr = 0.1,
                              sm = 'RR')

write_rds(event_rates_former, here("data_clean", "bayes_severity_former.rds"))

png(here::here('reports', 'figure', 'fig_8.png'), width=1024, height=546, res=120)
former_smoker_severity <- forest(event_rates_former,
                                        sortvar = Author,
                                        xlim = c(0.2, 10),
                                        rightlabs = c('RR', '95% CI', 'Weight'),
                                        leftlabs = c('Author', 'logRR', 'SE'),
                                        lab.e = 'Former smoker',
                                        lab.c = 'Never smoker',
                                        print.tau2 = F,
                                        col.diamond = 'blue',
                                        col.diamond.lines = 'black',
                                        col.square = 'black',
                                        col.square.lines = 'black',
                                        digits.sd = 2)
null <- dev.off()

event_rates_former

##              RR           95%-CI %W(random)
## Guan, Ni 2.9648 [1.7660; 4.9774]       13.0
## Rentsch  1.1283 [0.7885; 1.6146]       15.8
## Hadjadj  1.1111 [0.7419; 1.6640]       15.0
## Feuth    0.5312 [0.0702; 4.0199]        2.2
## Kuderer  1.6857 [1.3421; 2.1172]       17.9
## Gu       1.9805 [1.4715; 2.6657]       16.8
## Petrilli 1.2889 [1.1430; 1.4534]       19.1
## 
## Number of studies combined: k = 7
## 
##                          RR           95%-CI    z p-value
## Random effects model 1.5213 [1.1025; 2.0990] 2.55  0.0106
## 
## Quantifying heterogeneity:
##  tau^2 = 0.1372 [0.0115; 0.9881]; tau = 0.3704 [0.1074; 0.9941];
##  I^2 = 71.6% [38.4%; 86.9%]; H = 1.88 [1.27; 2.77]
## 
## Test of heterogeneity:
##      Q d.f. p-value
##  21.15    6  0.0017
## 
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - Sidik-Jonkman estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau

running_meta_count <- running_meta_count %>%
  full_join(., meta %>%
  rename("study" = author), by = "study") %>%
  select(study)

Figure 8: Former versus never smokers

0.5.5 Table 5: Smoking and COVID-19 mortality

#Table 5
table_5_word <- table_5 %>%
  mutate(., sample = sample_with_deaths) %>%
  mutate(., deaths_percentage = formatC(deaths/sample*100, digits = 2, format = "f")) %>%
  mutate(., death_current_smokers_percent = formatC(death_current_smokers/deaths*100, digits = 2, format = "f")) %>%
  mutate(., death_former_smokers_percent = formatC(death_former_smokers/deaths*100, digits = 2, format = "f")) %>%
  mutate(., death_current_former_smokers_percent = formatC(death_current_former_smokers/deaths*100, digits = 2, format = "f")) %>%
  mutate(., death_never_smokers_percent = formatC(death_never_smokers/deaths*100, digits = 2, format = "f")) %>%
  mutate(., death_never_unknown_smokers_percent = formatC(death_never_unknown_smokers/deaths*100, digits = 2, format = "f")) %>%
  mutate(., death_not_stated_percent = formatC(death_not_stated/deaths*100, digits = 2, format = "f")) %>%
  mutate(., recovered_percentage = formatC(recovered/sample*100, digits = 2, format = "f")) %>%
  mutate(., recovered_current_smokers_percent = formatC(recovered_current_smoking/recovered*100, digits = 2, format = "f")) %>%
  mutate(., recovered_former_smokers_percent = formatC(recovered_former_smoker/recovered*100, digits = 2, format = "f")) %>%
  mutate(., recovered_current_former_smokers_percent = formatC(recovered_current_former_smokers/recovered*100, digits = 2, format = "f")) %>%
  mutate(., recovered_never_smokers_percent = formatC(recovered_never_smoker/recovered*100, digits = 2, format = "f")) %>%
  mutate(., recovered_never_unknown_smokers_percent = formatC(recovered_never_unknown_smoker/recovered*100, digits = 2, format = "f")) %>%
  mutate(., recovered_not_stated_percent = formatC(recovered_not_stated/recovered*100, digits = 2, format = "f")) %>%
  select(lead_author, sample, recovered, recovered_percentage, recovered_current_smoking, recovered_current_smokers_percent,
         recovered_former_smoker, recovered_former_smokers_percent, recovered_current_former_smokers,
         recovered_current_former_smokers_percent, recovered_never_smoker, recovered_never_smokers_percent,
         recovered_never_unknown_smoker, recovered_never_unknown_smokers_percent, recovered_not_stated,
         recovered_not_stated_percent,
         deaths, deaths_percentage, death_current_smokers, death_current_smokers_percent,
         death_former_smokers, death_former_smokers_percent, death_current_former_smokers,
         death_current_former_smokers_percent, death_never_smokers, death_never_smokers_percent,
         death_never_unknown_smokers, death_never_unknown_smokers_percent, death_not_stated, death_not_stated_percent) %>%
  write_rds(here::here('data_clean', 'table_5_word.rds'))

The number of studies reporting on hospitalisation by smoking status is: 35

quality_table_5 <- table_5_word %>%
  left_join(., quality_rating, by = 'lead_author') %>%
  select(lead_author, overall_quality)

a <- table_5_word %>%
  mutate(Author = lead_author,
         Population = sample,
         recovered = paste(
           paste(recovered, 
                 as.integer(recovered_percentage), sep = " ("), "%)", sep = ""),
         recovered_current_smoker = paste(paste(recovered_current_smoking, 
                                                 (recovered_current_smokers_percent), sep = " (")
                                  , "%)", sep = ""),
         recovered_current_smoker = na_if(recovered_current_smoker, "NA ( NA%)"),
         recovered_former_smoker = paste(paste(recovered_former_smoker,
                                               (recovered_former_smokers_percent), sep = " (")
                                 , "%)", sep = ""),
         recovered_former_smoker = na_if(recovered_former_smoker, "NA ( NA%)"),
         recovered_current_former_smoker = paste(paste(recovered_current_former_smokers,
                                               (recovered_current_former_smokers_percent), sep = " (")
                                         , "%)", sep = ""),
         recovered_current_former_smoker = na_if(recovered_current_former_smoker, "NA ( NA%)"),
         recovered_never_smoker = paste(paste(recovered_never_smoker,
                                      (recovered_never_smokers_percent), sep = " (")
                                , "%)", sep = ""),
         recovered_never_smoker = na_if(recovered_never_smoker, "NA ( NA%)"),
         recovered_never_unknown_smoker = paste(paste(recovered_never_unknown_smoker,
                                              (recovered_never_unknown_smokers_percent), sep = " (")
                                        , "%)", sep = ""),
         recovered_never_unknown_smoker = na_if(recovered_never_unknown_smoker, "NA ( NA%)"),
         recovered_not_stated = paste(paste(recovered_not_stated,
                                    (recovered_not_stated_percent), sep = " (")
                              , "%)", sep = ""),
         recovered_not_stated = na_if(recovered_not_stated, "NA ( NA%)")) %>%
  mutate(deaths = paste(
           paste(deaths, 
                 as.integer(deaths_percentage), sep = " ("), "%)", sep = ""),
         death_current_smoker = paste(paste(death_current_smokers, 
                                                 (death_current_smokers_percent), sep = " (")
                                  , "%)", sep = ""),
         death_current_smoker = na_if(death_current_smoker, "NA ( NA%)"),
         death_former_smoker = paste(paste(death_former_smokers,
                                               (death_former_smokers_percent), sep = " (")
                                 , "%)", sep = ""),
         death_former_smoker = na_if(death_former_smoker, "NA ( NA%)"),
         death_current_former_smoker = paste(paste(death_current_former_smokers,
                                               (death_current_former_smokers_percent), sep = " (")
                                         , "%)", sep = ""),
         death_current_former_smoker = na_if(death_current_former_smoker, "NA ( NA%)"),
         death_never_smoker = paste(paste(death_never_smokers,
                                      (death_never_smokers_percent), sep = " (")
                                , "%)", sep = ""),
         death_never_smoker = na_if(death_never_smoker, "NA ( NA%)"),
         death_never_unknown_smoker = paste(paste(death_never_unknown_smokers,
                                              (death_never_unknown_smokers_percent), sep = " (")
                                        , "%)", sep = ""),
         death_never_unknown_smoker = na_if(death_never_unknown_smoker, "NA ( NA%)"),
         death_not_stated = paste(paste(death_not_stated,
                                    (death_not_stated_percent), sep = " (")
                              , "%)", sep = ""),
         death_not_stated = na_if(death_not_stated, "NA ( NA%)")) %>%
  select(Author, Population, recovered, recovered_current_smoker, recovered_former_smoker, recovered_current_former_smoker,
         recovered_never_smoker, recovered_never_unknown_smoker, recovered_not_stated, deaths, death_current_smoker,
         death_former_smoker, death_current_former_smoker, death_never_smoker, death_never_unknown_smoker,
         death_not_stated)

         
a$Author <- to_upper_camel_case(a$Author, sep_out = ", ")
a$Author <- plyr::mapvalues(a$Author,
                               from = cleaned_names,
                               to = correct_names)  

numeric_columns <- 'Population'

a <- flextable(a) %>%
    set_caption(caption = 'COVID-19 mortality by smoking status') %>%
  colformat_num(col_keys = numeric_columns, digits = 0, na_str = '-', big.mark = ',') %>%
  colformat_char(na_str = '-')

a <- set_header_labels(a,
                       Population = "Population with mortality",
                       recovered = "N (%)",
                       recovered_current_smoker = "Current smoker (%)",
                       recovered_former_smoker = "Former smoker (%)",
                       recovered_current_former_smoker = "Current/former smoker (%)",
                       recovered_never_unknown_smoker = "Never/unknown smoker (%)",
                       recovered_never_smoker = "Never smoker (%)",
                       recovered_not_stated = "Not stated (%)",
                       deaths = "N (%)", 
                       death_current_smoker = "Current smoker (%)", 
                       death_former_smoker = "Former smoker (%)", 
                       death_current_former_smoker = "Current/former smoker (%)",
                       death_never_unknown_smoker = "Never/unknown smoker (%)",
                       death_never_smoker = "Never smoker (%)",
                       death_not_stated = "Not stated (%)") %>%
  add_header_row(top = TRUE, values = c("","Recovered", "Died" ), colwidths = c(2, 7, 7)) %>%
  theme_booktabs() %>%
  fix_border_issues() %>%
  set_table_properties(width = 1, layout = 'autofit')
a

COVID-19 mortality by smoking status

		Recovered							Died
Author	Population with mortality	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Never/unknown smoker (%)	Not stated (%)	N (%)	Current smoker (%)	Former smoker (%)	Current/former smoker (%)	Never smoker (%)	Never/unknown smoker (%)	Not stated (%)
Chen	274	161 (58%)	5 (3.11%)	5 (3.11%)	-	-	-	151 (93.79%)	113 (41%)	7 (6.19%)	2 (1.77%)	-	-	-	104 (92.04%)
Zhou, Yu	191	137 (71%)	6 (4.38%)	-	-	-	-	131 (95.62%)	54 (28%)	5 (9.26%)	-	-	-	-	49 (90.74%)
Yang, Yu	52	20 (38%)	2 (10.00%)	-	-	-	18 (90.00%)	-	32 (61%)	-	-	-	-	32 (100.00%)	-
Borobia	2226	1766 (79%)	113 (6.40%)	-	-	-	-	1653 (93.60%)	460 (20%)	44 (9.57%)	-	-	-	-	416 (90.43%)
Giacomelli	233	185 (79%)	-	-	53 (28.65%)	132 (71.35%)	-	-	48 (20%)	-	-	17 (35.42%)	31 (64.58%)	-	0 (0.00%)
Yao	108	96 (88%)	1 (1.04%)	-	-	-	-	95 (98.96%)	12 (11%)	3 (25.00%)	-	-	-	-	9 (75.00%)
Carillo-Vega	9946	8983 (90%)	795 (8.85%)	-	-	-	-	8188 (91.15%)	963 (9%)	99 (10.28%)	-	-	-	-	864 (89.72%)
Heng	51	39 (76%)	6 (15.38%)	-	-	-	-	33 (84.62%)	12 (23%)	1 (8.33%)	-	-	-	-	11 (91.67%)
Chen, Jiang	135	NA (NA%)	-	-	-	-	-	-	31 (22%)	-	-	4 (12.90%)	-	-	27 (87.10%)
Heili-Frades	4712	4086 (86%)	210 (5.14%)	659 (16.13%)	-	-	3217 (78.73%)	-	626 (13%)	23 (3.67%)	161 (25.72%)	-	-	442 (70.61%)	-
Kim, Garg	2490	2070 (83%)	128 (6.18%)	481 (23.24%)	-	-	1461 (70.58%)	-	420 (16%)	22 (5.24%)	161 (38.33%)	-	-	236 (56.19%)	-
Al-Hindawi	31	15 (48%)	0 (0.00%)	10 (66.67%)	-	5 (33.33%)	-	-	16 (51%)	1 (6.25%)	12 (75.00%)	-	3 (18.75%)	-	-
Louis	22	16 (72%)	-	-	7 (43.75%)	-	-	9 (56.25%)	6 (27%)	-	-	3 (50.00%)	-	-	3 (50.00%)
Soto-Mota	400	200 (50%)	-	-	23 (11.50%)	-	-	177 (88.50%)	200 (50%)	-	-	25 (12.50%)	-	-	175 (87.50%)
Garibaldi	747	634 (84%)	36 (5.68%)	129 (20.35%)	-	-	-	469 (73.97%)	113 (15%)	6 (5.31%)	36 (31.86%)	-	-	-	71 (62.83%)
Docherty	13364	8199 (61%)	370 (4.51%)	1832 (22.34%)	-	4179 (50.97%)	-	1818 (22.17%)	5165 (38%)	214 (4.14%)	1350 (26.14%)	-	2105 (40.76%)	-	1496 (28.96%)
Kuderer	928	807 (86%)	38 (4.71%)	262 (32.47%)	-	425 (52.66%)	-	31 (3.84%)	121 (13%)	5 (4.13%)	64 (52.89%)	-	44 (36.36%)	-	2 (1.65%)
Ramlall	11116	10498 (94%)	-	-	2771 (26.40%)	7727 (73.60%)	-	-	618 (5%)	-	-	208 (33.66%)	410 (66.34%)	-	-
Wang, Oekelen	57	43 (75%)	-	-	14 (32.56%)	-	-	29 (67.44%)	14 (24%)	-	-	7 (50.00%)	-	-	7 (50.00%)
Martinez-Portilla	224	217 (96%)	-	-	7 (3.23%)	-	-	210 (96.77%)	7 (3%)	-	-	0 (0.00%)	-	-	7 (100.00%)
Cen	1007	964 (95%)	-	-	87 (9.02%)	-	-	877 (90.98%)	43 (4%)	-	-	1 (2.33%)	-	-	42 (97.67%)
Klang	3406	2270 (66%)	-	-	492 (21.67%)	-	-	1778 (78.33%)	1136 (33%)	-	-	301 (26.50%)	-	-	835 (73.50%)
Wang, Zhong	5510	4874 (88%)	247 (5.07%)	1083 (22.22%)	-	3544 (72.71%)	-	-	636 (11%)	28 (4.40%)	214 (33.65%)	-	394 (61.95%)	-	-
Miyara	338	211 (62%)	13 (6.16%)	58 (27.49%)	-	141 (66.82%)	-	-	46 (13%)	1 (2.17%)	23 (50.00%)	-	21 (45.65%)	-	-
Cepelowicz	255	209 (81%)	-	-	28 (13.40%)	181 (86.60%)	-	-	53 (20%)	-	-	18 (33.96%)	28 (52.83%)	-	-
Zeng	1031	866 (84%)	-	-	69 (7.97%)	-	-	797 (92.03%)	165 (16%)	-	-	36 (21.82%)	-	-	129 (78.18%)
Chen, Yu	1859	1651 (88%)	32 (1.94%)	54 (3.27%)	-	1565 (94.79%)	-	-	208 (11%)	13 (6.25%)	12 (5.77%)	-	183 (87.98%)	-	-
Garassino	190	124 (65%)	-	-	92 (74.19%)	32 (25.81%)	-	-	66 (34%)	-	61 (92.42%)	-	5 (7.58%)	-	-
Gu	884	864 (97%)	40 (4.63%)	250 (28.94%)	-	219 (25.35%)	-	-	20 (2%)	0 (0.00%)	14 (70.00%)	-	6 (30.00%)	-	-
Sigel	88	70 (79%)	-	-	37 (52.86%)	-	-	33 (47.14%)	18 (20%)	-	-	11 (61.11%)	-	-	7 (38.89%)
Nguyen	356	308 (86%)	-	-	91 (29.55%)	-	-	217 (70.45%)	45 (12%)	-	-	23 (51.11%)	-	-	22 (48.89%)
de Souza	8443	7826 (92%)	-	-	95 (1.21%)	-	7571 (96.74%)	160 (2.04%)	617 (7%)	-	-	47 (7.62%)	-	560 (90.76%)	10 (1.62%)
Mendy	532	663 (124%)	-	-	160 (24.13%)	-	-	502 (75.72%)	26 (4%)	-	-	10 (38.46%)	-	-	16 (61.54%)
Shi	256	210 (82%)	-	-	128 (60.95%)	-	-	82 (39.05%)	46 (17%)	-	-	26 (56.52%)	-	-	20 (43.48%)
Xie	619	591 (95%)	-	-	43 (7.28%)	-	-	548 (92.72%)	28 (4%)	-	-	8 (28.57%)	-	-	20 (71.43%)

save_as_docx(a, path = here('data_clean', 'Table_5.docx'))

0.5.5.1 Smoking and COVID-19 mortality meta-analysis

The studies included in meta-analysis are c(“al_hindawi”, “kuderer”, “miyara_medrxiv”, “chen_yu”, “gu”).

# Data --------------------------------------------------------------------
meta <- tibble('author' = table_5$lead_author,
               'recovered_current_smoker' = table_5$recovered_current_smoking,
               'recovered_never_smoker' = table_5$recovered_never_smoker,
               'death_current_smoker' = table_5$death_current_smokers,
               'death_never_smoker' = table_5$death_never_smokers,
               'recovered_former_smoker' = table_5$recovered_former_smoker,
               'death_former_smoker' = table_5$death_former_smokers) %>%
  filter(author %in% included_studies$lead_author)

meta$author <- to_upper_camel_case(meta$author, sep_out = ", ")
meta$author <- plyr::mapvalues(meta$author,
                               from = cleaned_names,
                               to = correct_names)  


# Current smoker mortality ------------------------------------------
event_rates_smoker <- meta %>%
  mutate(., Ee = death_current_smoker) %>%
  mutate(., Ne = (death_current_smoker+recovered_current_smoker)) %>%
  mutate(., Ec = death_never_smoker) %>%
  mutate(., Nc = (death_never_smoker+recovered_never_smoker)) %>%
  rename('Author' = author) %>%
  select(Author, Ee, Ne, Ec, Nc)

event_rates_smoker <- metabin(Ee,
                              Ne,
                              Ec,
                              Nc,
                              data = event_rates_smoker,
                              studlab = paste(Author),
                              comb.fixed = F,
                              comb.random = T,
                              method.tau = 'SJ',
                              hakn = F,
                              prediction = F,
                              incr = 0.1,
                              sm = 'RR')

write_rds(event_rates_smoker, here("data_clean", "bayes_mortality_current.rds"))

png(here::here('reports', 'figure', 'fig_9.png'), width=1024, height=546, res=120)
current_smoker_mortality <- forest(event_rates_smoker,
                                         sortvar = Author,
                                         xlim = c(0.1, 5),
                                         rightlabs = c('RR', '95% CI', 'Weight'),
                                         leftlabs = c('Author', 'logRR', 'SE'),
                                         lab.e = 'Current smoker',
                                         lab.c = 'Never smoker',
                                         print.tau2 = F,
                                         col.diamond = 'blue',
                                         col.diamond.lines = 'black',
                                         col.square = 'black',
                                         col.square.lines = 'black',
                                         digits.sd = 2)
null <- dev.off()

event_rates_smoker

##                RR            95%-CI %W(random)
## Al-Hindawi 2.6129 [1.0897;  6.2655]       26.3
## Kuderer    1.2394 [0.5189;  2.9604]       26.4
## Miyara     0.5510 [0.0800;  3.7976]       13.3
## Chen, Yu   2.7594 [1.7101;  4.4526]       32.1
## Gu         0.0920 [0.0002; 47.1721]        1.9
## 
## Number of studies combined: k = 5
## 
##                          RR           95%-CI    z p-value
## Random effects model 1.6652 [0.6836; 4.0565] 1.12  0.2616
## 
## Quantifying heterogeneity:
##  tau^2 = 0.5843 [0.0000; 10.7177]; tau = 0.7644 [0.0000; 3.2738];
##  I^2 = 29.2% [0.0%; 72.5%]; H = 1.19 [1.00; 1.91]
## 
## Test of heterogeneity:
##     Q d.f. p-value
##  5.65    4  0.2270
## 
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - Sidik-Jonkman estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
## - Continuity correction of 0.1 in studies with zero cell frequencies

Figure 9: Current versus never smokers

# Former smoker mortality ------------------------------------------
event_rates_former <- meta %>%
  mutate(., Ee = death_former_smoker) %>%
  mutate(., Ne = (death_former_smoker+recovered_former_smoker)) %>%
  mutate(., Ec = death_never_smoker) %>%
  mutate(., Nc = (death_never_smoker+recovered_never_smoker)) %>%
  rename('Author' = author) %>%
  select(Author, Ee, Ne, Ec, Nc)

event_rates_former <- metabin(Ee,
                              Ne,
                              Ec,
                              Nc,
                              data = event_rates_former,
                              studlab = paste(Author),
                              comb.fixed = F,
                              comb.random = T,
                              method.tau = 'SJ',
                              hakn = F,
                              prediction = F,
                              incr = 0.1,
                              sm = 'RR')

write_rds(event_rates_former, here("data_clean", "bayes_mortality_former.rds"))

png(here::here('reports', 'figure', 'fig_10.png'), width=1024, height=546, res=120)
former_smoker_mortality <- forest(event_rates_former,
                                         sortvar = Author,
                                         xlim = c(0.1, 5),
                                         rightlabs = c('RR', '95% CI', 'Weight'),
                                         leftlabs = c('Author', 'logRR', 'SE'),
                                         lab.e = 'Former smoker',
                                         lab.c = 'Never smoker',
                                         print.tau2 = F,
                                         col.diamond = 'blue',
                                         col.diamond.lines = 'black',
                                         col.square = 'black',
                                         col.square.lines = 'black',
                                         digits.sd = 2)
null <- dev.off()

event_rates_former

##                RR           95%-CI %W(random)
## Al-Hindawi 1.4545 [0.5500; 3.8467]        6.5
## Kuderer    2.0926 [1.4645; 2.9900]       44.0
## Miyara     2.1905 [1.2919; 3.7141]       21.3
## Chen, Yu   1.7367 [1.0224; 2.9500]       21.2
## Gu         1.9886 [0.7771; 5.0893]        7.0
## 
## Number of studies combined: k = 5
## 
##                          RR           95%-CI    z  p-value
## Random effects model 1.9767 [1.5393; 2.5384] 5.34 < 0.0001
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0038 [0.0000; 0.0796]; tau = 0.0618 [0.0000; 0.2821];
##  I^2 = 0.0% [0.0%; 2.6%]; H = 1.00 [1.00; 1.01]
## 
## Test of heterogeneity:
##     Q d.f. p-value
##  0.85    4  0.9310
## 
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - Sidik-Jonkman estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau

running_meta_count <- running_meta_count %>%
  full_join(., meta %>%
  rename("study" = author), by = "study") %>%
  select(study)

Figure 10: Former versus never smokers

0.5.6 Study quality assessment

#Table 6
fair <- nrow(filter(quality_rating, quality_rating$overall_quality == "fair" ))
poor <- nrow(filter(quality_rating, quality_rating$overall_quality == "poor" ))

The number of Fair quality studies are 32 with 140 Poor quality studies

The number of studies included in Meta-analysis are 26

0.6 Supplementary tables

0.6.1 Data on clinical or virological diagnosis

0.6.2 Data on studies stratifying smoking status